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CHAPTER  I 


INTRODUCTION 


To  put  the  problem  of  navigation  into  historical  perspective, 
let  us  quote  from  Singh's  "Great  Ideas  of  Modern  Mathematics:  Their 
Nature  and  Use."  ^ 

. . . trade  . . . grew  to  such  an  extent  that  it  transformed  the  whole 
life  of  Western  Europe.  It  is  true  that  the  development  of  trade 
led  to  a steady  growth  of  manufacture  as  well,  but  throughout  this 
period  (which  lasted  till  about  the  middle  of  the  eighteenth  century) 
trade  on  the  whole  dominated  manufacture . Thus  it  was  that  the  minds 
of  men  were  occupied  more  with  problems  connected  with  trade,  such 
as  the  evolution  of  safe  and  reliable  methods  of  navigation,  than 
with  those  of  manufacture.  Consequently,  while  the  two  ancient 
sources  of  power,  wind  and  water,  were  developed  at  a steadily 
accelerating  pace  to  increase  manufacture,  the  attention  of  most 
leading  scientists,  particularly  during  the  last  three  centuries  of 
this  phase,  was  directed  towards  the  solution  of  navigational  prob- 
lems. The  chief  and  most  difficult  of  these  was  that  of  finding 
the  longitude  of  a ship  at  sea.  It  was  imperative  that  a solution 
should  be  found  as  the  inability  to  determine  longitudes  led  to 
very  heavy  shipping  losses.  Newton  had  tackled  it,  although  with- 
out providing  a satisfactorily  practical  answer.  In  fact,  as  Hessen 
has  shown,  Newton's  masterpiece,  the  Principia,  was  in  part  an 
endeavour  to  deal  with  the  problems  of  gravity,  planetary  motions 
and  the  shape  and  size  of  the  earth,  in  order  to  meet  the  demands 
for  better  navigation.  It  was  shown  that  the  most  promising  method 
of  determining  longitude  from  observation  of  heavenly  bodies  was 
provided  by  the  moon.  The  theory  of  lunar  motion,  therefore,  began 
to  absorb  the  attention  of  an  increasing  number  of  distinguished 
mathematicians  of  England,  France,  Germany  and  America. 

If  we  replace  the  moon  by  a near  - earth  satellite,  optical 
observations  by  radio  signals  and  slide  rule  by  an  electronic  computer, 
we  have  today's  equivalent  navigational  problem.  Navigation  is  not 
restricted  to  just  ships  at  sea.  We  also  must  navigate  near  - earth 
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satellites,  planetary  and  lunar  proles  and  manned  space  capsules. 

Multivariate  regression  analysis  is  the  basic  mathematical 
method  of  solving  the  navigational  problem  from  electromagnetic  data 
corrupted  by  noise.  Regression  analysis  is  also  known  as  the  method 
of  least  squares  estimation  developed  by  Gauss.  Optimization  theory 
is  a collection  of  methods  to  aid  in  the  improvement  of  the  parameter 
estimation  of  a multivariate  regression  analysis.  The  improvement  is 
made  both  geometrically  and  analytically,  and  the  analytic  methods  will 
be  emphasized  here.  The  need  for  optimization  techniques  arises  from 
the  problem  of  correlation  among  the  parameters  in  a regression  analy- 
sis problem.  That  is,  an  error  in  one  parameter  will  compensate  for 
an  error  in  another.  Proper  attention  to  geometry,  such  as  the  proper 
positioning  of  measuring  instruments,  can  reduce  the  correlation  error. 
But,  if  general  analytic  optimization  methods  can  be  found,  much  data 
that  would  have  had  to  be  discarded, can  be  used,  and  hence  some  solutions 
can  be  saved.  Analytic  optimization  techniques  improve  the  accuracy,  in 
the  sense  of  reduced  biases,  of  the  parameter  estimation. 

The  next  chapter  on  Fundamental  Space  establishes  the  basic 
mathematical  structure  within  which  it  has  been  found  convenient  to 
study  multivariate  regression  analysis  problems.  Chapter  III  on  Opti- 
mization Theory  presents  the  basic  theorems  on  analytic  optimization 
techniques.  For  completeness,  two  regression  analysis  techniques  are 
also  given  in  Chapter  III.  The  last  chapter  investigates  the  feasi- 
bility of  the  methods  developed  in  Chapter  III  to  a particular  simu- 
lated problem  of  ship  navigation.  The  theoretical  methods  developed 
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in  Chapter  III  have  general  applicability  to  navigation  of  near  - earth 
satellites , planetary  probes  and  manned  space  flights,  and  as  a matter 
of  fact  to  general  multivariate  regression  analysis  problems.  A geome- 
tric optimization  result  is  also  presented  in  Chapter  IV.  Two  appendi- 
ces are  included,  the  first  on  generation  of  random  numbers  and  the 
other  on  simulation  of  satellite  orbits. 


CHAPTER  II 


FUNDAMENTAL  SPACE 

Quadratically  Normed  Finite  Dimensional  Linear  Space 
The  theory  of  optimization  to  be  developed  here  is  given  a 
firm  conceptual  foundation  by  rigorously  defining  a mathematical  struc- 
ture called  a Quadratically  Normed  Finite  Dimensional  Linear  Space , S, 
which  governs  the  collection  of  elements  called  vectors  or  points . The 
space  S is  then,  first,  a linear  space  satisfying  the  following  three 
axioms : 

1)  The  additive  structure  of  the  linear  space  is  an  abelian 

group.  That  is,  to  every  pair  of  vectors  x,  y,  there  exists 

a unique  vector  sum  such  that: 

(a)  If  x and  y are  in  S,  then  their  sum  x + y is  in  S 
( closure ) . 

(b ) For  all  x,  y,  z in  S,  x+(y  + z)  = (x+y)+z 
(associativity ) . 

(c)  There  exists  a unique  vector  zero,  9,  called  the  origin, 
such  that  x + 9 = Xi 

(d)  There  exists  a unique  vector  (-x)  such  that  x ■+■  (-x)  = 9, 
for  all  x. 

(e)  For  all  x,  y,  inS,  x + y = y-fx  (commutativity). 

2)  The  multiplicative  structure  of  the  linear  space  is  such  that 
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scalars  are  admitted,  as  operators  on  the  additive  group.  That 
is,  for  every  vector  x and  every  scalar  a from  a field  F,  there 
exists  a unique  vector  such  that: 

(a)  If  x is  in  S and  a is  in  F,  then  ax  is  in  S (scalar  multi- 
plication is  closed). 

(b)  For  all  a,  b in  F,  a(bx)  — (ab)x  (scalar  multiplication  is 
associative ) . 

(c)  lx  = x for  all  x,  where  1 is  the  multiplicative  identity 
element  from  the  field  of  scalars  F. 

3)  The  multiplicative  structure  is  related  to  the  additive  struc- 
ture by  the  distributive  laws.  That  is,  for  all  vectors  x,  y 
in  S and  all  scalars  a,  b in  F: 

(a)  a(x  4- y ) ^ ax  4-  ay  (scalar  multiplication  is  distributive 
with  respect  to  vector  addition). 

(b ) (a  ® b)x  = ax  4-bx,  where  a © b represents  the  sum  of  a 
and  b in  the  field  of  scalars  F (multiplication  is  dis- 
tributive with  respect  to  scalar  addition). 

In  our  work,  we  shall  be  concerned  only  with  the  field  of  sca- 
lars F being  the  real  numbers.  To  be  more  specific,  we  can  call  the 
space  a real  linear  space.  Whenever  the  term  "linear  space"  is  used,  we 
shall  mean  real  linear  space,  unless  otherwise  stated. 

A quadratic  norm  on  a linear  space  S is  a real  valued  function 
whose  value  at  x,  with  respect  to  a symmetric  positive  definite  matrix 
R is  denoted  by  j j xj | ^ , and  which  satisfies  the  following  four  axioms: 
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1) 

2) 

3) 


For  all  x in  S,  j|x 

INr  = ° lff  x = e 


For  all  x in  S and  all  a in  F 


* 


(positivity) . 

IHIs=  Nil  * llR 


(homogeneity) . 


k)  For  all  x,  y in  S,  ||x+y||R  =[|  x ||r  )|y  ||  R j, 

(triangle  inequality ) . 

A linear  space  S together  with  a quadratic  norm  on  S is  called 
a quadrat ically  normed  linear  space  iff  the  topology  is  defined  by  the 
distance  function 

d (x,  y)  = |jx-yl|  R , (l) 

vhere  d(x,  y)  satisfies  the  following  four  axioms  of  distance,  in  which 
case  S is  also  a metric  space  (regular  topological  space  with  a(A 
locally  finite  base).  The  real  valued  distance  function  d is  defined  on 
S x S such  that: 

\ 

l)  <i  (x,y)  = 0 

) (positivity) . 

2)  d (x,y ) — 0 iff  x = y 

3)  d (x,y)  =d(y,x)  (symmetry). 

*0  ^ (x,y)  = d(x,z  ) + d(z,y)  (triangle  inequality). 

We  must  show  that  the  distance  function,  as  defined  by  equation 
(l)  satisfies  the  four  distance  axioms,  l)  through  U).  Before  doing 
this,  we  shall  first  state  the  remaining  property  which  we  want  our  space 
to  have . 


The  space  S is  finite  dimensional,  iff  the  following  fundamental 


finiteness  axiom  is  satisfied: 
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There  exists  a finite  number  of  vectors,  v^,  vg,  . ..,  vn  in  the 
space  S such  that  every  vector  x in  S can  be  written  in  one,  and 
only  one,  way:  x — a.^  + a2Vg  «f  ...  + anvn. 

We  shall  define  a realization  of  the  quadratic  norm  over  a fi- 
nite dimensional  linear  space  as  follows : 

||xj|R  - (xTRx£,  (2) 

where  x is  a finite  dimensional  column  vector,  the  superscript  "T"  de- 
notes the  transpose,  the  indicated  square  root  is  positive,  R is  square 
of  the  same  dimension  as  x,  and  the  operation  between  the  vectors  and 
the  matrix  is  ordinary  matrix  multiplication. 

We  must  show  that  the  quadratic  norm,  defined  by  equation  (2) 
satisfies  the  four  axioms  for  a quadratic  norm,  l)  through  k).  For 
definiteness,  let  the  dimensionality  of  x be  n. 

Proof  l)  For  all  x in  S,  (x  Rx)  i 0,  since  by  definition,  a 
symmetric  matrix  R over  the  field  of  scalars  F is  positive  definite 
iff  the  quadratic  form,  defined  by  the  matrix,  is  also  positive  def- 
inite. A quadratic  form  is  positive  definite  if  its  values  are  posi- 
tive except  when  all  of  the  variables  are  zero.  The  quadratic  form  in 
x,  defined  by  the  matrix  is: 

T 

x Rx, 

T 

where  x — (x^_,  Xg,  ...,  x^)  . Thus  our  definition  for  the  quadratic 
norm,  which  is  the  positive  square  root  of  the  quadratic  form  in  x, 

satisfies  l). 

2)  If  x — 9,  all  components  of  the  vector  must  be  0 

T _ 

and  x Rx—  0. 
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Conversely,  if  xTRx  = 0 and  R is  positive  definite,  then  all  the 


variables  are  zero  and  x~Q.  Thus  2)  is  satisfied 
3)  II  ax 


R 


(ax)TR(axj2  -|a|(xTRx)2  ='a|  ]|xj jR, 
for  all  scalars  a in  F,  and  3)  is  satisfied. 

4)  We  begin  the  proof  of  axiom  4)  as  follows: 


x + y 


R 


(x  + y)  R(x  + y) 


J, 


- [(x  + y)T(Rx  + $r) 


- I (x  + y)  Rx  + (x  + y)  §y 


(xTRx  ) + (yTRx)  + (xTfy)  + (yT3/)j 


Nov, 

where 


x - Nig  + (y1®1)  + (*T$r)  +|[y;j  ® • (3) 

(h5  +ms>2  =iHr + 2 hi  R ik  ii  r+imIr  » " 

2 Hi  R|[y||  r = 2 (y\£- 

We  must  now  show 

(yTRx)  + (xTRy)  ^ 2(xTRxF(yT5y  )*• 

But  the  absolute  value  of  the  left-hand  side  becomes: 

j(yTRx)  + (xT5^)|  = l(yTRx)|  + |(xT§y)[  = 2 [(xTRy)[  , 

where  we  have  assumed  the  symmetry  property  (x^Ry)  — (y^Rx), 
which  may  be  easily  verified  by  expansion. 

Thus  we  must  show: 

|(xTRy)i  = (xTRx  )2(yTRyF. 

But  this  is  a generalized  version  of  the  Schwarz  Inequality  which  will 
be  proven  in  the  following 

Theorem  1 (The  Quadratic  Schwarz  Inequality). 

| (xTRy)[  xTRx)2(yTRy)i 
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The  equality  holds  if  x and  y are  linearly  dependent,  i.e.,  there 
exist  a-p  ag>  not  both  zero  in  F such  that  a^x  -f-  a^y  — 9. 


Proof . If  x or  y equals  9,  the  inequality  reduces  to  0=0. 
Assume  that  x and  y are  not  zero,  and  let  b be  a real  number. 

From  equation  (3)  we  have: 

II  x +•  by;!  p = |!x;|  p + ' (by  )TRx  ; + xTR(by)  4-  ||by|]  ^ 

= ||x||  | + 2b(xV)J+b^||y||^  0 . (4) 

The  =0  follows  from  axiom  l)  for  a quadratic  norm.  Thus  (1+ ) 
is  a non-negative  quadratic  equation  in  the  real  number  b.  In  order 
to  ensure  that  (4)  be  non-negative  for  real  b,  the  roots  will  have  to 
be  imaginary  or  equal.  But  this  implies  that  the  discriminant  of  a 


quadratic  equation  be  negative  or  zero.  Hence: 


or 


2(xT3jr) 
T 


4||: 


2 

R li 


1 y 


<x  sh|-!MI  H M 


■We  see  that  equality  holds  when  x =£  9 iff  the  discriminant  is 
zero.  The  discriminant  is  zero  iff  the  roots,  which  must  be  imaginary 
or  equal,  are  equal  since  a zero  discriminant  gives  no  chance  for  an 
imaginary  root  to  exist.  Thus  the  equality  will  hold  iff  there  is  a 
real  value  of  b for  which  equation  (4)  will  equal  zero.  That  is,  for 
some  b,  1 x + byj!  = 0 iff  x + by  - 9 by  axiom  2)  for  a quadratic  norm. 
Let  b = a^/a^  where  a-^  £ 0 and  a2  are  real  numbers.  Then  x -f  (a^/a^)y 
= 9 implies  a^x  + agy  — 9,  which  implies  x and  y are  linearly  dependent. 
Also,  if  x — 9,  the  equality  in  the  theorem  holds  and  x and  y are  linear- 


ly dependent , 
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The  completion  of  the  proof  of  the  Quadratic  Schwarz  Inequality 
also  completes  the  proof  that  the  quadratic  norm,  defined  in  equation 
(2),  satisfies  the  four  axioms  for  a quadratic  norm.  Restating  axiom 
4),  in  terms  of  our  defined  quadratic  norm,  we  have  proven  the  follow- 
ing: 


Theorem  2. 


(x  -f  y)  R(x  -f  y) 
which  will  also  be  called  the  Quadratic  Minkowski  Inequality. 


2 = (xTRx)^  + (yTRy  )^, 


It  remains  to  show  that  the  distance  function  defined  in 
equat ion  ( 1 ) , 

d(x,y)  =||x  - y||  = (x  - y)TR(x  - y)  2 

_ L 

satisfies  the  four  distance  axioms  l)  through  4). 


Proof . From  axiom  l)  for  a quadratic  norm: 

H*-r||HS  0. 

r=» 

Hence,  d(x,y)  = 0,  and  the  distance  axiom  l)  is  satisfied.  From  axiom 
2)  for  a quadratic  norm: 

||  x - y ||  R = 0 iff  x - y = 9. 

Hence,  d(x,y)  — 0 iff  x — y,  and  the  distance  axiom  2)  is  satisfied. 

From  axiom  3)  for  a quadratic  norm: 

||  X - y[|  Il(-i)(y  - x)j|  r=  f-lf  II  y - x||R  = |[y  - x|]  R. 
Hence,  d(x,y)  — d(y,x),  and  the  distance  axiom  3)  is  satisfied.  From 
axiom  4)  for  a quadratic  norm  we  have: 

lix  " y||  R = 11  x - z - y zll  R =llx  - zli  R +IU  - yll  R- 

Hence,  d(x,y)  = d(x,z)  -f-  d(z,y),  and  the  distance  axiom  4)  is  satisfied. | 
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In  some  cases  we  shall  be  interested  in  a quadratic  inner  pro- 
duct space.  Thus,  a real  linear  space  S with  vector  elements  x,y  is 
a quadratic  inner  product  space  if  there  is  defined,  on  SxS,  a real 
valued  function  (x,y)^,  called  the  inner  product  of  x and  y with  re- 
spect to  R,  which  is  a positive  definite  symmetric  matrix  with  elements 
from  the  field  F of  real  scalars,  satisfying  the  five  axioms: 

!)  (x+y>z)R  = (x,z  )R  -f  (y,z  )R,  (linearity), 

2)  (x,y)R=  (symmetry), 

3)  (ax,y)  - a(x,y)R,  a from  the  real  scalar  field  F 

(homogeneity). 


k)  (x,x)  = 0 
R 

5)  (x,x)  = 0 iff  X = 9 

A 


(positivity ) . 


We  shall  define  the  quadratic  inner  product  in  terms  of  the 
quadratic  norm.  Thus  if  S is  a quadratically  norned  linear  space  with 
vector  elements  x,  y,  we  define  the  inner  product  with  respect  to  the 
symmetric  positive  definite  matrix  R on  SxS  as  a real  valued  function 


(5) 


With  this  definition,  letting  y — x,  we  can  replace  the  last  two  axioms 
by  the  single  axiom 


**')  (*>*)  "'ll  x!l  u • 

R " 11  R 

Using  our  definition  for  a realization  of  the  quadratic  norm  in 


equation  (2), 


R 

let  us  apply  this  to  the  definition  for  the  quadratic  inner  product  in 
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equation  (5), 


(x,y)  - £ 
R 


x + yll  R -||  x - y|| 


and  obtain  a realization  for  the  quadratic  inner  product . 
Thus  the  right-hand  side  becomes: 


1 

k 


(x  +y)T  R (x  +y)  - (x  - y)T  R (x  - y) 


- £ (xTRx)  + (yTRx)  + (x\)  + (yT$r)  - (xTRx)  + (yTRx)  + (x^)  - (yT^), 
= £ [ 2(yTRx)  + 2(xT^/  = £ [U(xTRy)^  . 

Thus  a realization  of  the  quadratic  inner  product  of  any  two  vectors  x, 
y,  with  respect  to  the  symmetric  positive  definite  matrix  R,  correspond- 
ing to  our  defined  quadratic  norm,  equation  (2)  is: 

(x,y)  = (xTRy).  (6) 

n 

We  note  that  if  y = x,  we  have 


(x,x)  - (xTRx)  — 1|  x 


R 


R 


so  that  our  quadratic  inner  product  correctly  reduces  to  axiom  b')  for 
the  quadratic  inner  product  space  S.  With  definition  (6)  for  the  quad- 
ratic inner  product,  we  must  show  that  the  remaining  quadratic  inner 
product  axioms  are  satisfied, 

(x  +y,z)R  = [(x  -f  y)V  = (xTRz)  + (yTRz)  = (x,z)R  +-  (y,z)R, 


satisfying  l).  Then, 


(x>y)R-  [xTRy]  - | yTRx'  - (y,x)  , 


I 

which  satisfies  2). 

(&x,y)R  = r(ax)TRyJ  = a(xTRy)  = a(x,y) 

_ J R 

for  all  a from  the  scalar  field  F.  Thus  3)  is  satisfied  and 

(xT§y) 
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is  a realization  of  a quadratic  inner  product  of  x and  y with  respect 
to  the  symmetric  positive  definite  matrix  R. 

The  following  theorem  now  may  be  proved  easily.  It  may  be  vis- 
ualized in  the  Euclidean  plane  as  a relation  between  the  major  and 
minor  diagonals  of  the  parallelogram  with  its  sides,  and  is  called  the 


Quadratic  Parallelogram  Theorem.  For  all  elements  x,  y,  in  a 
quadratic  inner  product  space  S,  defined  by  the  quadratic  norm  in 
equation  (5),  the  following  holds: 


X + H!  R+ll*  - HI  R ~ 2 H X 


2 _ 


R 


+ 2||y 


Proof . By  axiom  4')>  the  left  side  above  becomes 
(x  + y,  x -t-y)R  + (x  - y,  x - y)R  = (x,  x + y )R  +-  (y,  x +•  y)R  +■  (x,x  - y)R 

+ (■  y>  x > y) 

K 

- (x  3-y,x)R  + (x  fy,y)R  +■  ( x - y,x)R  - (x  - y,y)R 

= (x>x)R  + (y>x)R  + (x>y)R  + (y>y)R  +■  (x>x)R  +■  (-y>x)R  ' (x>y)R 

- (-  y,y) 

- 2(x,x)r  -+  2(y,y)R  - 2 j|  x ||  R - 2 ||  y ||  R 

by  repeated  application  of  axioms  l),  2),  3)  and  1+'). 


Quadratic  Finite  Dimensional  Pre -Hilbert  Space 
A quadratically  normed  finite  dimensional  real  linear  space 
will  be  called  a Quadratic  Finite  Dimensional  Pre -Hilbert  Space  iff 
the  quadratic  norm  satisfies  the  additional  condition  given  in  the 
Quadratic  Parallelogram  Theorem: 

ilx  "HI  R +llx  - HI  R - 2(  INI  R + IMI  R ) 
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Since  our  definitions  for  a realization  of  the  quadratic  norm 
(2)  and  of  the  quadratic  inner  product  (6)  satisfy  the  theorem,  our 
space  S is  a Quadratic  Finite  Dimensional  Pre-Hilbert  Space. 

Quadratic  Orthogonality 

We  shall  restate,  for  convenience,  the  realization  of  the  quad- 
ratic inner  product  of  any  two  vectors  x,y,  with  respect  to  the  symmet- 
ric positive  definite  matrix  R,  corresponding  to  our  defined  quadratic 
norm: 

(x,y)R  - (*TRy). 

The  following  identities  then  result  from  the  definition: 

(x,y)R  = (x,$r)  = (Rx,y)  - xTRy . 

We  now  may  make  the  following 


Definition  1.  Suppose  x and  y are  vectors  in  a quadratic  fi- 
nite dimensional  pre-Hilbert  space  S and  suppose  R is  a symmetric  pos- 
itive definite  matrix  with  elements  from  the  field  F of  scalars  . 

Then  the  vectors  x and  y are  orthogonal  vectors  with  respect  to  B.  or, 

in  short,  quadratic  orthogonal  vectors,  iff  (xty)„  = Q.  We  may  symbol- 

R 

ize  this  by  writing  x J_  y. 

R 

From  this  definition,  we  see  that  when  R — I,  the  identity  ma- 

T 

trix,  we  obtain  the  ordinary  definition  for  orthogonality,  (x  y)  — (x,y) 

= 0.  When  R is  a diagonal  matrix  D of  rank  r,  say  D — diag  (d  ,d_...,d  ), 

1 2 r 

T 111 

then  it  is  possible  to  write  D - P P where  P = diag  (d^2,  d^2 . ,d  2). 

(x,y)D  = (xT5y)  = (xTPTPy)  =0. 


Thus 
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Now  let  u — Px  and  v — Py.  Then  the  above  becomes: 

(x,y)  - (u,v)  - (uTv)  =:  (u,v)  - 0. 

- I 

It  can  be  show:?  that  for  every  symmetric  positive  definite  matrix  R, 

T 

there  exists  a nonsingular  matrix  P such  that  R - P P,  so  that  the 
above  relations  hold  for  all  admissible  R matrices. 

We  next  wish  to  define  the  quadratic  orthogonality  of  two  sub- 
spaces, for  which  we  in  turn  require  the  definitions  for  linear  sub- 
space, dimension,  linear  basis,  and  linear  dependence. 

Definition  2.  The  vectors  x^,  x.  f . . . in  a linear  space  X 

are  linearly  dependent  vectors  iff  there  exists,  from  the  field  F of 
scalars,  a^,  (i  — l,...,m),  not  all  zero  in  F,  such  that  a^x^  -(- e^Xg 
■+-•••+  amxn  0*  These  vectors  are  linearly  independent  iff  they  are 
not  linearly  dependent. 

Definition  3-  A linear  basis  (or  a coordinate  system)  in  a 
linear  space  X is  a set  V of  linearly  independent  vectors  in  X such 
that  every  vector  in  X can  be  written  as  a linear  combination  of  vectors 
in  V. 


Definition  4.  The  dimension  of  a finite  dimensional  linear 
space  X is  the  number  of  elements  in  a linear  basis  of  X.  (This  defi- 
nition follows  from  the  theorem  that  the  number  of  elements  in  any  lin- 
ear basis  of  a finite  dimensional  linear  space  is  unique.) 


Definition  5 • A non-empty  subset  W of  a linear  space  V is  a 


1 6 


linear  subspace  (linear  manifold)  iff  for  every  pair,  x and  y,  of  vec- 
tors in  W,  every  linear  combination  ax  + by,  vhere  a,b  are  from  the 
field  of  scalars,  is  also  contained  in  W.  Thus  the  linear  manifold  is 
itself  a linear  space. 

We  are  now  ready  to  define  vhat  we  mean  by  quadratic  orthogo- 
nal subspaces . 

Definition  6.  Suppose  S and  S are  linear  manifolds  of  a 

X J 

quadratic  finite  dimensional  pre-Hilbert  space  S and  suppose  R is  a 
symmetric  positive  definite  matrix.  Then  the  linear  manifolds  Sx  and 
Sy  are  orthogonal  linear  manifolds  with  respect  to  R or,  quadratic  or- 
thogonal linear  manifolds,  iff  (x.y)  =0  for  all  x in  S and  v in  S . 

R X * y 

We  may  symbolize  this  by  writing 

It  is  easy  to  define  a quadratic  orthogonal  set  of  vectors. 

Definition  7 • Suppose  < |>  is  a set  of  vectors  in  a quadrat- 

ic finite  dimensional  pre -Hilbert  space  S and  suppose  R is  a symmet- 
ric positive  definite  matrix.  Then  the  setjx^  is  an  orthogonal  set 
of  vectors  with  respect  to  R or,  quadratic  orthogonal  set  of  vectors, 
iff  every  pair  of  vectors  in  < x^  | is  orthogonal  with  respect  to  R, 
that  is,  (xfXj)R  = 0 for  all  i,  j such  that  i i j. 

It  is  easy  to  prove  the  following  generalization  where  S and  R 
are  defined  as  above: 


Quadratic  Pythagoras1  Theorem.  If  x and  y are  vectors  in  S 
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and  R is  a symmetric  positive  definite  matrix,  such  that  x | y, 


then 


T,  -1- 


vhere  we  define  [|z||  = (zTRz)2  for  all  z in  S. 

R 


Proof. 


x +y  ' r “ Sx  + (x  + y) 

- (xTRx)  + (yTRx)  + (xT§y)  + (yTRy) 


L 

= 1 1 x 


p + (y*x)  + (*>y)w  +!|y||  \ 

K R R 11  11  R 


where  we  used  the  realization  for  the  quadratic  inner  product.  From 

the  condition  in  the  theorem,  x | y implies  (x,y)  — 0,  and  by  symmetry 

_ R 


R 

of  the  quadratic  inner  product,  (y,x)  - 0.  Thus: 

R 

2 ~ii  n 2 
R 


x+yii  * =[MI  l + h\\ 


I 


Generalized  Quadratic  Pythagoras'  Theorem.  If  the  set  <|  x 


of  vectors  are  in  S,  and  R is  a symmetric  positive  definite  matrix 
such  that 


'■i_[_Rxj>  for  aH  i>  j - l,...,n,  and  i i j, 


then 


i — 1 


n p n 

T.  r --  I 


R 


i — 1 


Proof . From  the  condition  in  the  theorem,  (x^x.)  — 0,  for  all 
i,  j - 1, . . .,n  and  i t j.  | 

In  the  following  definitions,  suppose j x^  is  a set  of  vectors 
in  a quadratic  finite  dimensional  pre-Hilbert  space  S and  suppose  R 
a symmetric  positive  definite  matrix. 


is 
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Definition  8.  The  set  j is  an  orthonormal  set  of  vectors 
respect  to  fl  or,  quadratic  orthonormal  set  of  vectors  iff  i) 
every  pair  of  vectors  in  j xij>  is  orthogonal  with  respect  to  R,  and  ii) 
each  vector  is  normal  with  respect  to  R,  that  is 

!|  xi||  R - (xi>xi)R  - (xi^i)  = 1,  for  all  i. 

Now,  every  linear  basis  c(x.l  of  a linear  manifold,  S of  S can 

l XJ  a 

be  orthonormalized . That  is,  we  can  find  a set  jv^j  In  terms  of  the 
linear  basis  jx^j>  of  linearly  independent  vectors  in  Sa  such  that  eve- 
ry vector  in  Sa  can  be  written  as  a linear  combination  of  vectors  in 
K}  and  such  that  |vi|  are  orthonormal  with  respect  to  R.  We  shall 
state  this  as  a theorem  which  will  use  a proof  that  leads  to  the  Quad- 
ratic Gram-Schmidt  Orthonormalizing  process. 

Theorem_3.  Supposejx^ j>is  a finite  linear  basis  consisting  of 
k elements  in  a linear  manifold  S&  of  a quadratic  finite  dimension- 
al pre -Hilbert  space  S,  and  suppose  R is  a symmetric  positive  defi- 
nite matrix.  Then,  we  can  find  constants 

all 

a21  a22 

a.,  a a 

31  32  33 


such  that  the  elements 
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V1  - allxi 


v2  ~ a21xl  + a22x2 
v - a x 4-  a x„  + a x 

3 31  l 32  2 33  3 


are  orthonormal  with  respect  to  R. 

Proof.  We  normalize,  vith  respect  to  R,  the  first  vector  by- 
dividing  it  by  its  quadratic  norm: 

V1  = X1  i V1  -V!lvlll  R- 

Now  define  a vector  w^  orthogonal  (with  respect  to  R)  to  v^  by  set- 


ting 


Wo  XP  “ (xOt  vi  ) 


1 /RV1 


then  normalize. 


V — v 
2 2 


w. 


2"  R 


Next  define  a vector  w orthogonal  (with  respect  to  R)  to  v and  v. 


Wo  ~ x - (x  ,v  ) v - (x  ,v  ) 


v 


— w / w 


3 3 ' 3'  2'R  2 ' 3'  l'R  1 3 3'  11  "3"  R 

The  general  form  for  the  jth  vector,  j = k,  is  given  as  follows: 


Wj  Xj  ‘ 


J - 1 


i - 1 


(xj^vi)Rvi  J VJ  " Vllwj 


R 


It  is  thus  seen  that  each  w^ , and  hence  each  v^  is  a linear  combination 


of  the  | x^j>,  and  since  the|x^  > is  a set  of  linearly  independent  vectors, 
each  w.  i 0 and  hence  the  quadratic  norm  j|w  ||  ^ 0 by  positivity  of 

quadratic  norms . We  must  show  that  the  vectors  are  normal,  with 
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respect  to  R: 


(,Wr  ' 


Vi 


WJ  ' 


IWJ||  R inn  a 


R 


T „ 

uj  ?WJ 


VJ| 


"jllR 


w 


= 1, 


j — 1 , • » • ,k . 

To  complete  the  proof  ve  shall  show  that  w.  is  orthogonal  (with  respect 

J 

to  R)  to  v]_,  v2,  hy  induction.  Thus,  as  an  induction  basis, 

show  (w2,  v]_)^  ~ 0 • That  is 

\ \ — * ( \ _ l ( \r  t r ir  tr  \ 


( [x2  - (x2»  vi)rV1J  ' vl)R  ' <x2>  V1>R  ‘ l(x2»vl  Vl^vl)_ 

~(.x2>  vi)r  ' (x2'  vlVVl»  V1}R 

= (V  VA  ' (X2>  V1  }R  ° ' 

since  (v1,  v^  = 1. 

Assume  the  induction  hypothesis  (w  , v , ) — 0 for  all  m^j  and  i^m. 

m i R 

Then  for  i ■=  j , 

("o’  vi>R  = 


/ t1 

\ \ 

h - £ 

q-l 

(Y  VrT4)  ’ Vi) 

J-l 

z (xr  Vr  ‘ ^ °V  Vs(V  vi5r 

q-1 


- (xj’  Vr  - (xJ’  Vr  - 0 ’ 

since  (w^,  ) - 0 implies  (vm,  vi)R  - 0 for  i“=m,  and  (vi,  — l.| 

We  now  need  to  define  the  concepts  of  a dimensionally  complete 
orthonormal  set  of  vectors  and  an  orthonormal  linear  basis,  where  we 
recall  the  previous  definitions  of  S and  R above  . 


Definition  9.  The  set  <j  xi  j>  of  vectors,  orthonormal  with 
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respect  to  R,  is  dimensionally  complete  in  the  space  S iff  the  only  vec- 
tor y in  S which  is  quadratically  orthogonal  to  every  vector  of  j x.  > 
is  the  vector  y — 9.  A dimensionally  complete  set  of  vectors  which  is 
also  orthonormal  with  respect  to  R,  is  called  a dimensionally  complete 
orthonormal  set  with  respect  to  g,  or,  a d- complete  quadratic  orthonor- 
mal  set . 


Definition  10.  A set  ij  x^  j>  of  vectors  in  a linear  manifold  S& 
of  S is  a quadratic  orthonormal  linear  basis  iff  the  set  is : 

1)  A quadratic  orthonormal  set  of  vectors  (with  respect  to  R). 

2)  d-complete  in  the  linear  manifold  Sa . 


Quadratic  Orthogonal  Projection 

Definition  11.  Suppose  T is  a subset  (not  necessarily  a linear 

manifold)  of  a quadratic  finite  dimensional  pre-Hilbert  space  S,  and 

suppose  R is  a symmetric  positive  definite  matrix.  Then  a vector  s 

in  S and  the  subset  T are  orthogonal  with  respect  to  R,  or,  quadratic- 

ally  orthogonal,  iff  (s,t)  — 0,  for  all  t in  T.  We  may  symbolize  this 

R 

by  writing  s | ^T,  or  by  s | ^t , for  all  t in  T.  The  set  of  all  such  vec- 
tors s is  the  quadratic  annihilator  (with  respect  to  R)  of  T,  and  will 
be  denoted  by: 

T — — |s  in  S : (s,t)  — 0 for  all  t in  T >. 

Note  that  s | T implies  and  is  implied  by  s in  T 
We  also  make  the  following  definitions: 


-hR 
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(T  *L^"L§)  r T J§i 


Theorem  U . If  T is  any  subset  of  a quadratic  finite  dimen- 
sional pre-Hilbert  space  S and  R is  a symmetric  positive  definite 

J-R 

matrix,  then  the  quadratic  annihilator  T - is  a linear  manifold  of 

S. 


Proof.  From  the  definition  of  linear  manifold,  the  subset 


-4* 


-41 


— is  a linear  manifold  iff  for  every  pair  s^,  s 2 of  vectors  in 


T -,  every  linear  combination  as1  + bs2,  a,  b from  the  field  of  sea- 


f R 

lars,  is  also  contained  in  T From  the  definition  of  quadratic  an- 


nihilator, asx  + bs2  is  contained  in  T if 

(aB1  + bs2,  t)R  = 0 

for  all  t in  T . From  the  axioms  for  a quadratic  inner  product  space  we 
have : 

(as]_  + bs2,  t)R  = (aSl,  t)R  + (bs2,  t)R 

“ a(sx,  t)R  + b(e2,  t)R 

- aO  -f  bO  - 0,  for  all  t in  T, 

-*-R 


since  s^  and  s2  are  members  of  T — 


I 


Theorem  5 . If  T and  U are  subsets  of  a quadratic  finite  dimen- 
sional pre -Hilbert  space  S and  R is  a symmetric  positive  definite 
matrix,  then: 

1)  TCT  ^ 

2 ) If  TC  U,  then  U^CT^ 
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3) 


Proof:  To  facilitate  the  proof,  we  first  observe  the  follow- 
ing implications  from  the  definition:  If  X and  Y are  subsets  of  S 

then 

a)  x in  X implies  x i ^X  — , since  x in  X implies  x | Ry  for 


all  y in  X — 


h)  y | rY  implies  and  is  implied  by  y in  Y "J"- 
Thus  to  prove  l),  suppose  t in  T.  This  implies  t | T by 
).  Ify  b),  the  last  implies  t in  (T  — ) — . Hence  t in  T implies  t 


a 


in  T^^,  that  is,  TCT  - -. 

To  prove  2),  x in  U ~ implies  from  condition  b)  above,  x i U. 

-l-R 

Hence,  since  x | and,  if  TC  U,  then  also  x j T implies  from  condition 
b)  above,  x in  T Hence  x in  U _J"-  implies  x in  T that  is,  U "i"- 

CT  -,  provided  that  TCU. 

For  the  proof  of  3)  we  observe  from  l),  TC T and  by  2) 

this  implies  (T  — ) — C T — . Also  from  l),  replace  T by  T — on 

both  sides,  giving  T — c These  two  inclusions  imply 

-*-R  — *JR  -i-R  -j-R  i p . -p  i p 

T - T ~ . We  made  use  of  the  identity  (T  - — 

- (T  5)  by  applying  the  quadratic  annihilator  to  (T  ~-) 

- T - - giving  (T  -)  - - = T £ h £ = (T  - -)  -,  where  the 


last  equality  comes  from  the  definition. 
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Definition  12.  Suppose  T and  U are  subsets  of  a quadratic  fi- 
nite dimensional  pre-  Hilbert  space  S and  R is  a symmetric  positive 
definite  matrix^  then  the  subsets  T and  U are  quadratic  orthogonal 


2k 


subsets  (with  respect  to  R ) iff  (t,u)  =0  for  all  t in  T and  all  u 

R 

in  U,  and  may  be  symbolized  by  T , U.  From  this,  it  readily  follows 
that  T , pT  and  | pT. 

Decomposition  Theorem.  If  T is  a complete  linear  manifold  of 
a complete  quadratic  finite  dimensional  pre -Hilbert  space  S,  and  R 
is  a symmetric  positive  definite  matrix,  then  the  quadratic  anni- 
hilator  of  T,  T — , is  called  the  quadratic  orthogonal  complement 
of  T.  Any  vector  x in  S can  be  decomposed  uniquely  in  the  form: 
x = t + r,  where  t in  T and  r in  T 

Proof.  To  prove  uniqueness,  assume  that  the  decomposition  is 
not  unique: 

x “■  ^1  ri  " tg  -f  r2>  with  t-^  and  tg  in  T and 

_ in 

rx  and  rg  in  T -,  t tg  and  r±t  rg. 

We  can  form  the  new  vector 

v ~ t2  ~ tl  ~ rl  “ r2*  vhere  (t2  " ti^  in  T and  (ri  “ r2)  in  T "A'“* 
From  T_^T  - we  see  that  w t Rw  implies  (w,w)  - 0 implies  w = 9. 

Hence  we  conclude  that  tg  — t^  and  r2  — and  the  decomposition  is  u- 
nique . 

WTe  must  now  show  the  existence  of  the  decomposition.  Let|v.|> 
be  a quadratic  orthonormal  linear  basis  in  T.  Let  t a^  where 
ai  ~ (x>vi)R*  x = t 4-  r.  This  equation  defines  r,  so  that  all 

that  remains  is  to  show  that  r is  quadratically  orthogonal  to  T.  That 
is,  we  must  show  that  (r,v^)^  ~ 0 for  all  i.  Solve  the  defining  equa- 
tion above  for  r . Then  from  the  axioms  for  a quadratic  inner  product 
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space  we  have : 


(r,Vj)R  - ((x  -Z  .jVjl.v  ) - (x,Vj)s  - (IVl>T  ) 


JR 


- (*,^)R  - (Z  (x.v^v^v,) 


R i J R 


- ^'Vr  • (x'vl)„(vl.'vj)B  - (x'v2)b(v2'v.1) 


R' 


J V V a V • J • • • 

R 2'  J'R 


(x>v )„  - ••• 

J r J Jr 


— (x,v^ *1—0,  for  all  j where  use  is  made  of 


/ _ r _ Jl,  1 = j 

the  identity  (v.,v  ) - c)  -\0,  i t j.  Thus  we  have  shown  that 

A J R !j 

(r,v.)  - 0 for  all  j and  therefore  r is  quadratically  orthogonal  to  T.| 

J £_ 


The  Decomposition  Theorem  leads  to  the  following  definitions  of 
orthogonal  projections: 


Definition  13 . Suppose  T is  a complete  linear  manifold  of  a 
complete  quadratic  finite  dimensional  pre-Hilbert  space  S,  and  R is  a 
symmetric  positive  definite  matrix.  Given  x in  S and  suppose  x — t -f  r 
where  t in  T and  r in  T1^.  The  vector  t is  defined  as  the  quadratic 

orthogonal  projection  of  x on  T (with  respect  to  R)  and  the  vector  r 

-LR 

Is  defined  as  the  quadratic  orthogonal  projection  of  x on  T — , the 
quadratic  orthogonal  complement  of  T. 

This  leads  to  the  definition  of  the  projection  operator: 

Definition  Ik . Suppose  T is  a complete  linear  manifold  of  a 
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complete  quadratic  finite  dimensional  pre -Hilbert  space  S and  R is  a 
symmetric  positive  definite  matrix.  For  each  x in  S we  denote  the 
quadratic  orthogonal  projection  operator  (vith  respect  to  R)  onto  T as 
PT(R),  and  the  quadratic  orthogonal  projection  onto  T of  x will  be  de- 
fined by 


PT(R)X  - t,  where  t in  T. 

We  may  obtain  an  analytic  definition  for  the  projection  oper- 
ator in  terms  of  a linear  basis,  as  demonstrated  in  the  following  two 
theorems . 


Theorem  6.  Let  T be  a complete  one -dimensional  linear  mani- 
fold of  a complete  quadratic  finite  dimensional  pre-Hilbert  space  S 
and  let  R be  a symmetric  positive  definite  matrix.  For  each  x in  S 

let  P . be  the  quadratic  orthogonal  projection  onerator  (with  re- 
T(R) 

spect  to  R)  onto  T and  let  s be  a non-zero  vector  in  T.  Then: 


PT(R)X  ■ 8'  vhere  8 ln  T' 

’ R 


Proof . Select  a scalar  b so  that  bs  is  the  quadratic  ortho- 
gonal projection  of  x on  T . That  is,  from  simple  geometric  consider- 
ations, we  have 

bs_^(x  - bs). 

This  implies 


Expanding, 


(bs,  (x 

the  above  becomes 

(bs,x) 


- bs)) 
R 

+ (bs, 


= 0. 


- bs)  = 0 
R 


1 
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2 

b(s,x)  - b (s,s)  = 0. 

K K 

Thus  b = (x's)g  , 

(6'6)r 

(X's)R  ■ 

giving  P / -.x  - bs  - _ — _=  s where  s in  T.  5; 

(s,s)R  B 

We  see  immediately  from  the  above  theorem  that  if  the  vector  s 
is  normalized  vith  respect  to  R then  the  above  projection  becomes: 
^T(r)x  — ® in  T,  (s^s)^  — 1. 

Notice  that  when  R = I,  the  identity  matrix,  the  above  projection  re- 
duces to  ordinary  orthogonal  projection. 

From  the  fundamental  finiteness  axiom,  we  know  that  there  ex- 
ists a finite  set  of  basis  vectors js^j  such  that  any  vector  in  our 
space  can  be  written  as  a linear  combination  of  the  basis  vectors. 

This  fact,  combined  with  the  last  theorem,  proves  the  following 


Quadratic  Orthogonal  Expansion  Theorem.  Let  the  set  of  non- 
zero vectors  j j>  be  a quadratic  orthonormal  linear  basis  of  a 
linear  manifold  T of  a complete  quadratic  finite  dimensional  pre- 
Hilbert  space  S and  let  R be  a symmetric  positive  definite  matrix. 
For  each  vector  s^  let  T^  be  the  corresponding  complete  one  dimen- 
sional linear  manifold  of  T such  that  s in  T^  For  each  arbitrary 
vector  t in  T let  .(r)  Ot^d-TS-tic  orthogonal  projection  oper- 

ator (with  respect  to  R ) onto  T^ . Then  t may  be  represented  as 


follows : 
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^ £ ^ip  (r)^  £ ®i  in  Tp  • 

i 1 - i 

This  theorem  leads  to  the  following 

Definition  15 . Suppose  all  the  definitions  in  the  above  theo- 
rem apply  here . Then 

£ (t,s.)R  s. 
i 

is  called  the  Quadratic  Fourier  Series  for  t,  and  the  constants 


are  the  Quadratic  Fourier  Coefficients  of  t. 

For  completeness,  ve  shall  define  the  projection  operator  onto 
T1^,  below. 

Definition  15.  Suppose  T is  a complete  linear  manifold  of  a 
complete  quadratic  finite  dimensional  pre -Hilbert  space  S and  R is  a 
symmetric  positive  definite  matrix.  The  quadratic  orthogonal  pro- 
lg.ction  operator  onto  T_^§  is  denoted  by  , the  quadratic  ortho- 

gonal complement  of  T.  The  quadratic  orthogonal  projection  onto  T"^, 
of  an  arbitrary  vector  x in  S will  be  defined  by 

pT-lR  x ~ rj>  where  r in  T^. 

The  following  theorem  contains  some  of  the  properties  of  the 
quadratic  orthogonal  projection  operator. 

Theorem  7 (Quadratic  Orthogonal  Projection  Operator). 

Let  T be  a complete  linear  manifold  of  a complete  quadratic  finite 
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dimensional  pre-Hilbert  space  S and  let  R be  a symmetric  positive 

definite  matrix.  For  each  x in  S let  P be  the  quadratic  ortho- 

T(R  j 

gonal  projection  operator  (with  respect  to  R)  and  let  x = t + r 
where  t in  T and  r in  T1^.  Then  has  the  following  properties: 

x)  pT(R)t  ~ b for  a11  t in  T, 

2)  PT(R)r  ~ 9 f°r  a11  r in 

^PT(R)X,y^R  ^ ^X,PT(R)y^R"  (syTnmetry  Property  of  pT^R^), 

11 ) PT(R)(PT(R)X)  “ Pt(r)x>  (indempotent  property  of  PT(R)), 

PT(R)^X  + ~ PT(R)X  + PT(R)y * for  a11  x>  y in  S> 

PT(R)(bx)  ~ bPT(R)x>  vhere  b is  a scalar  from  the  field. 


7) 


pt(r)xII  r = 


R ’ (Pt(r)  1s  bounded)* 


Properties  5)  and  6)  indicate  Prp^R^  is  a linear  operator, 


Proof,  l)  and  2)  are  immediate  consequences  of  the  definition. 

To  prove  3),  ve  use  the  identities  x = P . ^x  + PT^  and  y 

= PT(R)y  n°te  that  (PT(R)X)VPT-Ry)  implleS  (PT(R)X'PT^y)R 

“ 0 and  similarly  with  (P^)  ^(P  y). 

(PT(R)X^y)R  * (PT(R)X'  PT(R)y  4’PT-#)R 

(PT(R)X,PT(R)y )r  + (PT(R)X>  P^  y )r  = (PT(R)X'PT(R)y ^ 

(pT(R)x,  pt(r))r  + (p  [plR  x'  PT(R)y)R 

(pT(R)x  + PT-^X>  PT(R)y^R  ~ PT(R)y^R 


To  prove  k),  we  have,  from  the  definition,  that  P , x = t. 

T(R) 

Thus  we  must  show  that  the  left-hand  side,  P , N(pm/  .x)  = P , s(t) 

* T(R)V  T(r)x'  rT(R)^'> 
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becomes  equal  to  t.  But  this  follows  from  property  l)  above. 

To  prove  5),  suppose  x - t -f  r,  y — s + p,  where  t,  s in  T and 
. hR 

r,  p m T . Then  x + y — (t  + s)-f(r+p)  vhere  (t  -f-  s)  in  T and 
(r  + p)  in  T^.  Thus 

pT(?)(x  + y)  - t + s = pt(5)X  + PT(g)y. 

To  prove  6)  multiply  x - t -f  r by  the  scalar  b to  give 
bx  - bt  -f  br.  Then  Pij(R)(bx)  - bt  - bPT^R^(x)  where  we  used  the  def- 
inition pT(R)U)=t. 

Proof  of  7):  For  all  x 

HX'i  R ^lPT(R)X  + PTi*  X'l  R 

^PT(R)X  + PTJ?X'  PT(R)X  + Pt"LPx ^R 
^PT(r)X*  PT(R)X^R  + PT'^pX^R 

~ llPT(R)xll  R II  PTJ^XH  R > 

IHU  = ^pt(r)xIIr  (7) 

Define  ||PT(g)||  = ^nf  ^ it? , vhere  u = |m:  ||PT(g)x||  * g M ||x||  X (8) 

Thus,  from  (7)  and  (8)  we  have : 

llPT(R)H  = 1-  | 

% analogy , the  quadratic  orthogonal  projection  operator  of  x 
onto  T~ii,  Pt-lR,  can  be  shown  to  have  equivalent  properties,  l)  through 
7). 


xneorem  8.  Let  T,  S,  R,  be  defined  as  in  the  Quadratic 

Orthogonal  Projection  Operator  Theorem;  and  let< s . | T p . , bp 

( 1 y T (R)  De 
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defined  as  in  the  Quadratic  Orthogonal  Expansion  Theorem.  For  each 
x in  S let  x = t + r where  t in  T and  r in  T^.  Then 

PT(R)X  - t 1 PT1(R)X  - E (x>si)R  si  • 

PT(S)X  = PT(R)^‘  + = PT(B)t  + PT(Rf  = ‘ + 8 = t frOT> 

l)  and  2)  of  the  above  theorem.  Also, 

^ PTi(R)x  ~ E PTi(R)^t  +r)  ‘ E PT.(R)t  E PTj  (R)r 

l 1 i 1 ~ i 1 ~ 

’ ^ PTi(K)t’ 

“d  E (x.^)R=i-  E (t  + r.sp^r  £ (t,st)Rsi+  £ (r,.^ 

i - i i ~ i 

= E (‘-“iVi  • 

i 

The  conclusion  in  the  theorem  follows  by  the  result  in  the 
Quadratic  Orthogonal  Expansion  Theorem: 

' ‘ £ PT1(R)t  ' E V 81  in  Ti-  I 

In  order  to  prove  the  next  theorem,  ve  must  be  assured  of  the 
existence  of  a vector. 

Definition  17 . Suppose  S&  is  a complete  subset  of  a complete 

quadratic  finite  dimensional  pre -Hilbert  space  S.  A subset  S of  S is 

a 

convex  iff  whenever  x and  y are  in  SQ,  so  is  the  line  segment  that 
joins  them:  mx  -+  (l  - m)y  in  SQ,  0 = m = 1. 

We  see  that,  if  Sa  is  a complete  linear  manifold,  it  is  convex, 
since  by  definition  ax  + by  in  Sa  for  all  x,  y in  and  all  a,  b in 
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the  scalar  field. 

The  important  theorem  of  this  section  now  follows  and  is  called: 

Quadratic  Orthogonal  Projection  Theorem. 

Let  T be  a complete  linear  manifold  of  a complete  quadratic  finite 
dimensional  pre-Hilbert  space  S and  let  R be  a symmetric  positive 
definite  matrix.  Let  the  set  of  non-zero  vectors  <|s.\  be  a quad- 
ratic orthonormal  linear  basis  of  T and  let  be  the  quadratic 

orthogonal  projection  operator  (with  respect  to  R)  onto  T.  Let  x 
be  an  arbitrary  vector  in  S.  There  exists  a yQ  in  T such  that 


x - y. 


R 


= - y 


R 


for  any  y in  T.  The  minimum  is  attained  iff 


T(R) 


x - yr 


Proof.  From  the  Decomposition  Theorem,  any  vector  x in  S can 
be  decomposed  uniquely  in  the  form 

i P 

x — t + r,  where  t in  T and  r in  T — . 

From  the  last  theorem  we  saw  that 

PT(R)X  ~ t ' 

Now,  for  any  y in  T,  j|x  - y||  2 - (x  - y,x  - y) 

° R 


- ! [x  - t]  + [t  - y]  , X - tj  + rt  - y] 


But  x - t - r,  (x  - t ) J j-.  T,  where  t,  y and  (t  - y) 
Hence  the  foregoing  becomes: 


7* 


are 


- (x  - t,  x - t)  + (t  - y,  t - y)  , 
■K  R 


in  T. 
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since  (t  - y,  x - t)  - 0 because  (t  - y ) | R (x  - t).  The  above  nay 


be  written 


x - yll  r - llx  • -til  l +11*  - yll  r- 


Since  each  quantity  is  positive. 


x - t 


II  * s H= 


2 

R 


|x  - y 

demonstrates  the  existence  of  the  minimum  value  for  the  quadratic 
norm.  We  may  write  (9)  above  in  the  following  form,  where  we  have 
made  the  substitution 

t = E pTl(R)*  = E (t^L  Bi>  si  in  TiCi  T, 


(9) 


i i 

and  y — E where  the  b^  are  arbitrary  scalars,  where 

i 

quadratic  orthonormal  linear  basis  of  T: 


is  a 


x - E (t,3l)R  6l]|  * § || x - EblBi 


2 

R' 


Let  us  expand  the  right-hand  side: 
2 


biSil-  R 


(x  ~ E bisi>  x " E bisi) 


i i a 

— (x  Rx)  - (b^s^,x)^  - E (xjb^s^)  + E (bisi>bjs j )r 

i i - y 

-!IxI!r  - Ebi(si^)R  - E bi(x»si)R+  E bibj  (8i>6j) 

i i y 


R 


-II xil  r - 2 Ebi(si>xE  + Eb?  +E  (si>x)R  -E  (6i»x)R 


- IIX||  R + E 0>i  -(Si,x)R)2  - E (8i>X) 


- i i - i 

2 v / v2 
R 


(10) 


i i 

Thus  we  see  that  the  minimum  is  attained  iff 

b±  = (8i>x)R  (11) 

which  we  defined  previously  as  the  quadratic  Fourier  coefficients  of  x, 


3^ 


since 


R 


Z_i  (s^ >x  )R  > 


l!x  - ? b±sill  I s II  >■ 

1 i 

all  quantities  are  positive,  and  equality  is  attained  if  condition  (ll) 

is  satisfied. 


R 


Since  x-t  + r where  t in  T and  r in  T^,  we  have 

bi  = (si^)p  = (s.,  t + r)^  = (8i,  t)R  + (Si,r)R  = (S±,t) 
since  s^  in  T,  r in  T"1?  implies  s^  | ^ implies  (s^,r)  — 0. 

Multiply  through  by  s.  and  sum: 

y = \ Vi  = 3 (Wj*!  - £ %l(5,t  - t = pi(5)x, 

the  last  equality  coming  from  the  previous  theorem.  Thus  we  have  dem- 
onstrated that  the  minimum  is  attained  iff 


T(R) 


x - y- 


I 


Corollary  (Quadratic  Bessel's  Inequality). 

__,(si>x)p  = j|x;|Rp  where  / s1  ^ is  an  orthonormal  linear 


R ~ II  “ Hr 

basis  with  respect  to  R 


Proof, 


inf 
b . 


^ x ' t'blSl11  l ■ V(Sl'X,R' 

which  is  obtained  from  equation  (10)  in  the  theorem  setting  b 
Then,  since  all  quantities  are  positive,  solving  for 

Zj  (si>x)p 

i ~ 

gives 


- (s,,x) 


R 


s ||,|||. 
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For  purposes  that  will  appear  later,  we  want  to  place  the  form 
on  the  left-hand  side  of  equation  (10)  (in  the  theorem  proof)  into  ma- 
trix form.  First,  suppose  that  x is  an  nXl  vector,  and  let  the  sum- 
mation extend  from  1 to  q.  Thus  is  the  i^*1  component  of  a qxl 
vector  and  is  the  i^1  column  vector  of  an  nxq  matrix,  where  is 

an  nXl  vector.  Thus  we  can  write: 


x - 


V 
L 

i = 1 


biSi  = 


fxl 

1! 

X2! 

. > 


nJ 


*1 

b2 


> 


N 


/ 

X1 

S11  s12  • • • slqj 

bl 

X2 

S21  6 22  s2q; 

b2 

~ < 

• 

• 

• 

lxn. 

> " 

. : 

• • • 

• • • 

snl  sn2  snq 

' > 

Kl 

— x - Sb 

where  S represents  the  nxq  matrix  and  b the  qxl  vector. 


CHAPTER  III 


OPTIMIZATION  THEORY 

Introduction  to  Optimization 

Optimization  theory  is  concerned  with  the  determination  of  the 
precision  of  the  best  estimate  of  a set  of  parameters , perhaps  by  simu- 
lation studies,  and  the  utilization  of  this  information  to  improve  the 
parameter  estimation  of  a multivariate  regression  analysis.  The  pre- 
cision of  the  best  estimate  of  the  set  of  parameters  is  defined  statis- 
tically in  terms  of  a variance -covariance  matrix.  If  our  multivariate 
regression  analysis  results  in  a minimum  variance  estimate  it  is  pos- 
sible to  formulate  a distinction  between  estimation  theory  and  opti- 
mization theory:  In  estimation  theory  a minimum  variance  estimate  is 

achieved  by  properly  incorporating  the  covariance  matrix  of  the  obser- 
vation errors;  in  optimization  theory  an  improved  estimate  is  obtained 
by  proper  utilization  of  the  covariance  matrix  of  the  parameter  errors. 

In  optimization  theory  we  are  not  greatly  concerned  with  the 
distribution  function  of  the  parameters  being  estimated  because  what- 
ever the  distribution  is,  we  transform  it  into  a uniform  distribution. 3 
Thus  the  covariance  matrix  of  the  parameters  relates  to  the  ellipsoid 
of  concentration  (or  error)  and  hence  represents  the  precision  of  the 
estimate  of  the  parameters.  This  measure  of  precision  is  then  utilized 
in  two  ways  to  improve  the  parameter  estimation:  l)  Geometric 
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optimization,  and  2)  analytic  optimization. 

As  an  example  of  geometric  optimization,  suppose  that  we  want 
to  determine  the  optimum  array  of  instruments  to  use  for  determining 
space  position  of  a particular  trajectory  or  orbit  of  a satellite.  By 
calculating  the  error  ellipsoids  for  various  combinations  of  instru- 
ments at  specific  points  along  the  trajectory,  it  is  possible  to  deter- 
mine optimum  arrays  for  each  point.  However,  even  with  modest  numbers 
of  instruments  the  problem  can  quickly  get  out  of  hand.  For  example, 
if  we  were  given  15  instruments  and  wanted  to  determine  the  error  el- 
lipsoids for  all  possible  combinations  of  two  up  through  eight  instru- 
ments at  a time  for  one  trajectory  point,  the  total  number  of  ellip- 
soid calculations  would  exceed  22,000.  Even  with  high  speed  computers 
this  load  is  intolerable,  and  so  a method^  has  been  developed  to  re- 
duce the  number  of  computations  to  about  1000 per  point  for  this  prob- 
lem. The  method  selects  up  to  the  20  best  solutions  at  each  stage  of 
k instruments  and  calculates  the  ellipsoids  at  the  k + 1 stage  by 
adding  one  new  instrument  at  a time.  The  best  solutions  are  deter- 
mined by  a function  of  the  error  ellipsoids.  A related  method  was  de- 
veloped to  determine  the  optimum  instrument  positions  for  an  entire 
trajectory  or  satellite  orbit.  Such  optimization  problems  are  studied 
by  simulation  of  orbits  and  instruments. 

In  order  to  study  analytic  optimization  techniques  we  must 
first  turn  our  attention  to  the  regression  analysis  techniques  of  pa- 
rameter estimation.  In  regression  analysis  problems  it  is  not  uncommon 
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for  all  the  parameters  to  become  correlated  with  each  other.  That  is, 
an  error  in  one  parameter  will  compensate  for  an  error  in  another,  and 
the  estimates  will  be  said  to  be  biased.  Four  methods  have  been  de- 
rived to  improve  the  biased  estimates,  and  the  related  theory  is  pre- 
sented in  this  chapter.  Each  method  has  been  tested  for  its  effective- 
ness, and  the  results  on  a simulated  navigation  problem  are  given  in 
the  following  chapter.  The  four  methods  are  called:  l)  Uncorrelated 

Estimation,  2)  Quadratic  Pseudoinverse  Estimation,  3)  Constrained  Esti- 
mation and  4)  Functionally  Dependent  Constrained  Estimation. 

Two  regression  analysis  techniques  are  presented:  l)  Multiva- 

riate Regression  Analysis  and  2)  Sequential  Regression  Analysis.  The 
first  technique  is  also  called  a weighted  least  squares  estimation. 
Sequential  estimation  is  related  to  communication  engineering,  and  is 
derived  directly  from  the  first  technique . 

In  all  our  estimation  procedures  we  assume  a linearized  model. 

If  the  equations  governing  the  relations  between  measurements  and  param- 
eters to  be  estimated  are  not  linear,  then  we  assume  that  a Taylor 
series  expansion  up  to  linear  terms  is  valid.  We  also  assume  that  the 
errors  associated  with  the  measurements  from  an  instrument  are  uncor- 
related random  variables  forming  a normal  Gaussian  distribution.  We 
assume  all  matrices  that  are  to  be  inverted  can  be  inverted  to  what- 
ever accuracy  we  desire. 

This  chapter  concludes  with  three  special  sections  which  are 
needed  for  the  previous  sections.  They  are  titled  Derivation  of  Co- 
variance,  Probability  Ellipsoid  of  Error  and  Matrix  Calculus. 
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Multivariate  Regression  Analysis 5 
In  the  following  we  shall  suppose  that  the  space  S governing 
the  vectors  is  a complete  quadratic  finite  dimensional  pre -Hilbert 
space  and  that  R is  a symmetric  positive  definite  matrix.  The  reali- 
zation of  the  quadratic  norm  of  a vector  x (with  respect  to  R)  is  de- 
fined by 

]'i  x ||  R = (xTRx)2. 


Suppose  that  y is  a linear  function  of  q unknown  parameters 

"til. 

bx,  b2,  b , so  that  for  the  i value  of  y we  can  write  the  exact 

expression 


?i  * ailtl  + ai2b2  + * • • + aiqbq'  vhere  the  coefficients 
aij  are  known.  Suppose  that  a sequence  of  measurements  y^,  ...  , yn 
has  been  made,  and  define  the  measurement  vector  as 

yT  =jyi,  y2<  •••  » yn}- 

lxn  1 J 

Define  the  vector  of  parameters  by 


and  let  A 

n>rq 


T _ \ 

b — s b-^,  bg,  . . . 
1 X q 1 

with  rank  q,  q<cn. 


) 


Theorem  1.  The  parameter  estimate  b for  the  absolute  minimum 
of  the  quadratic  minimizing  function, 


Q(b)  - y - Ab 


-II  R> 

where  R is  called  the  weighting  matrix,  is  given  by 


nxn 


^ T -1 
b - (A  RA)  A Ry. 


4o 


Proof.  Apply  the  differential  operator  matrix  V.  to  the  quad- 
ratic  norm  function  C(b),  set  equal  to  zero,  evaluate  for  b - b and 

A 

solve  for  b.  Thus  the  condition  is 

First  expand  by  using  the  definition  for  the  realization  of  the  quad- 
ratic norm: 

Vb  [(y  - Ab)TR  (y  - Ab)]  = 0, 

Vb  [yTPy  - yTRAb  - bTATRy  - bTATRAbj  - 0. 

Then,  since  y and  R are  independent  of  b and  R is  symmetric, 
-[(V/)PAb  + ( VbbTATR)y]  - [(VbbTATR)y  + (VbyT)RAb] 

+ [(VbbTATR)Ab  + (\7bbTAT)RAb]  = 0. 

Then,  carry  out  the  operation  and  evaluate  for  b = b,  where  VbyT  - 0: 

- 2ATRy  + 2ATRAb  - 0 . 

Therefore,  b - (ATRA)"1ATR  y . (l) 


V.QCS)  = Vjy  - Ab| 


2 

R 


We  must  verify  that  this  value  of  the  estimated  parameters 

actually  gives  an  absolute  minimum.  Nov,  by  definition,  the  coefficient 

of  y in  equation  (l)  is  the  quadratic  useudoinverse  0 of  A with  re- 

n 

spect  to  R , 
nxn 

A*(H)  = (ATRA)_1ATR  > 

qxn  ~ 

where  it  can  be  shown  that 


(a*r),*r)  ' 

Thus  we  may  write  (l)  above  in  the  form 

= AfR)y  , 


See  Section  on  Quadratic  Pseudoinverse. 


(2) 


or  equivalently 


y-Ab, 

since  substituting  (2)  into  (3)  gives 

y=Ab-AA*y  - y,  where  AA  — I 

— l.~ i\  J — r — n ’ 

- - nxn 

or  substituting  (3)  into  ( 2 ), 

Ab 


(3) 


8 = 4 


A ^ 

- b , where  A A — I 

“S’  ql 


To  show  that  (l)  is  an  absolute  minimum,  first  expand  the  quadratic 
minimizing  function  Q(b): 


Q(b)  - ||  y - Ab 


R 


PO 


rn  T T T T 

- y Ry  - y RAb  - bA^+bA  RAb . 

Now  substitute  (l)  in  the  equivalent  form  of  equation  (3), 

y = A&, 

into  (4 ) above : 

C(b ) ~||  yj|  R “ ^TATRAb  - bTAT RAS  + bTATRAb  4 bTATRAS  - bTATRA$ 


“II  yll  I 4 (b  - b)TATRA(b  - b)  - bTATRAb 


Q(b)  -!|  y|j  ^ 4 1|  A(b  - b)||  ^ - ]j  AS 


— 1 IR  ‘ 


(5) 


Thus  we  see  from  expression  (5)  that  Q(b)  attains  a minimum  value  iff 

(6) 


A 

b - b 


as  given  by  (l),  since  from  (5), 


y - Ab 


- Il  R 


= , y 


1 2 

R 


a,  p 

-bn  R » 


all  quantities  are  positive,  and  equality  is  achieved  when  (6)  is 
satisfied. 

The  covariance  matrix  corresponding  to  the  estimate  for  the 

/S 

parameters,  b,  is  given  in  the  following 


k2 


Theorem  2.  Let  S represent  the  covariance  matrix  correspond- 

v 

ing  to  measurement  errors,  and  let  S ^ represent  the  covariance  ma- 

b 

trix  corresponding  to  parameter  estimate  errors.  Suppose  the  param- 
eter estimate  b is  given  by 

b - (ATRA)“1AT^ 

with  all  matrices  and  vectors  defined  as  in  the  previous  theorem. 
Then  >p 


T -1  T T -1 

Sc  — (A  RA)  A RS  RA(A  RA)  . 

y 


(7) 


Proof . Defining  as  before  the  quadratic  pseudoinverse  of  A 
with  respect  to  R as 

A*  - (ATRA)_1ATR, 

• R 

allows  us  to  write  the  parameter  estimate  in  the  form 

A.  * 

b - A y . 

-R 

From  the  section  on  Derivation  of  Covariance,  the  covariance  matrix  of 

A 

b is  given  by 

M.  rp 

(8) 


q x q q n X n n >rq 


5 6 = sy 


Now, 


*>T 


„T 


-l.T  i T 


.T, 


i 


*-i 


(A^)  = (A  RA)  Oi  ~ ?A(A  RA)  , 


- „T 


where,  since  R is  symmetric,  R — R , and  where  use  is  made  of  the  fact 


that  the  transpose  of  an  inverse  is  the  inverse  of  the  transpose.  Sub- 
stituting the  above  result  and  the  definition  for  A*  into  (8)  yields 
the  result  identical  to  (j)  in  the  theorem.  | 


The  following  corollary  readily  follows: 


^3 


Corollary  1 . Let  the  weighting  matrix  R be  chosen  such  that 


S'"  - §y. 

Then 

b - (a  Sy 

and 

It 

- -y 


Weighting  the  quadratic  minimizing  function  by  S , results  in  what 

v 

is  called  in  the  literature  a minimum  variance,  or  Markov,  esti- 
mate . 


Proof . Substitute  Sy  for  R in  equation  (7)  in  the  above 

theorem. 

g 

Sequential  Regression  Analysis0 
Sequential  regression  analysis  requires  the  same  foundation 
as  that  for  the  main  theorem  in  the  section  on  Multivariate  Regression 
Analysis.  In  this  type  of  analysis,  the  measurements  are  processed 

A 

one  at  a time  as  they  are  received.  The  parameter  estimate  b is  then 
updated  as  each  measurement  is  processed.  This  type  of  estimation  is 
particularly  well  suited  for  space  navigation  in  real  time  since  it  is 
not  necessary  to  invert  large  sized  matrices. 

Suppose  that  we  have  obtained  the  estimate  of  the  parameters 
b from  a sequence  of  n measurements.  From  the  main  theorem  in  the 
Multivariate  Regression  Analysis  section  we  know  that  the  solution  is 
given  by 

-N  rp  -1  rn 

b - (A  RA)  A~R  y , 


kk 


vhere  the  quadratic  minimizing  function  is 

Q(b)  - ] y - Ab  , . 

When  we  add  a new  measurement  to  Q(b)  we  shall  write 

Q(b)  -||  y - Abi!  £ 

where  R is  now  of  rank  n 4 1, 

7 ~ ; yi#  y2*  ” ' > v yn  _ 2} 

and  A is  (n  4-  l)  xq,  of  rank  q.  The  solution  is  then  given  by 

Cl  nw  . — . -]  rp 

b ~ (d  3 A)  A1  R y . 

If  we  define  R,  where  r is  a scalar,  as 

r 1 

n - R 0 u , 

_ ' > which  is  usual  for  a minimum  vari- 

|_0  r J 

ance  estimate , then  the  solution  becomes 

where  a - j an  + 1^1,  an  + 2 , • • • •>  an  + p q } are  known  constants . 

This  form  requires  the  inversion  of  a (qxq)  - matrix.  To  avoid  this 
we  require  first  two  lemmas. 


A 

b 


- , A RA  4 


1-1 


a ra 


A R y 4 a r y 


n 4 


(9) 


Lemma  1. 


aPa^  4 r“b 


■i  _ r T t 1 -i  t 
a - I 14  a raP  x a ra. 


where  P'1  - ATR  A 


Proof.  Multiply  both  sides  by  I 4 a raP  , giving 


T "j  t r *p 
I 4 a raP  | a , aP  a 4 r 


-1  1 -1  _ T 

1 a — a ra 


T _1  T rp 

a r r 4 a1  raP  a1 


[ aP  a^  4 r"b 


-1  - T 

a - a ra 


T -1  T1 

a r ! r x 4 aP  a1 


aP  a*  + P-1  : -la  = aT 


ra . 
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Matrix  Inversion  Lemma.  Suppose  P — 


P'1  - ATRA.  Then  P = P - PaT 


-1  T 1 -1 
P -f  a ra  where 


T -1 
aP  a1  + r 


-1 


aP. 


Proof.  Apply  Lemma  1 to  the  right  hand  side  of  the  above 
equation,  giving 

P*  = P - P 


T 

I 4 a raP 


-1  T _ 
a raP  . 


Now,  P can  be  written  | P 1 ! 1 and  using  the  knowledge  that  the  product 


of  inverses  is  the  inverse  of  their  products  gives 

P = P 


-1  T 
P + a ra 


-1  T 
a raP. 


Multiply  both  sides  by 


-1  T 
P + a ra 


to  give 


— 1^ 

P P - 


-1  T 
P •+•  a ra 


P - a raP 


_ _ -1  t T 

I — P P + a raP  - a raP  . 


I 


The  following  theorem  then  states  the  form  of  equation  (9)  for 
which  we  seek. 


Sequential  Estimation  Theorem.  Suppose  an  estimate  of  the  pa- 
rameters after  n measurements  is  given  by 
b - (ATRA)_1  ATR  y , 
and  after  n -4  1 measurements, 


A, 

b-  = 


T T 

A RA  4 a ra 


-1 


T T 

A R y 4 a r y 


n 4-  1 


Suppose  P_1  - ATRA  and  P = ATRA  4 aT  ra . Then 


1 

T 

A 

aP 

a r 

y , - ab 

J 

L n 4 1 J 

Proof.  We  begin  the  proof  by  writing  the  equation  in  the  normal 


form, 


T T 

A RA  + a ra 


~ T T 

b - A Ry  + a r y 


n + l' 


T 

A RA 


A m 

b — AiR  y . 


Subtract  the  second  equation  from  the  first  to  obtain 

~ a r t 1~_t 

(b  - b)  -f  ara  b-a  r ) . . , 

n -+-  1 


T 

A RA 


~ A A ^ A a 

Define  b - b — db  , so  that  b — db  -f-  b,  and  make  both  of  these  substi- 


tutions into  the  equation  giving 


r .t  . i a 

' T 

/\ 

T 

A RA  j db  + 

a ra 

db  4- 

a ra 

0 - T 
b — a r y 


n + 1 


Solve  for  db : 


A 

db  - 


T T 

A RA  + a ra 


n>  m a 

ary  . - a rab 
n -f  1 


The  quantity  in  the  left  hand  bracket  is  P;  apply  the  Matrix  Inversion 

£ A 

Lemma  and  re substitute  b - b for  db  to  give 
* = $ + ! £ - PaT(aP  aT  + r-1)’1^  | aTr  [ yn  + ± - ab]  . | 


Un correlated  Estimation 

Optimization  theory  uses  the  precision  in  the  estimation  of  the 
parameters  of  a multivariate  regression  analysis  to  improve  the  esti- 
mation. The  multivariate  regression  will  be  restricted  to  that  which 
results  in  a minimum  variance  (Markov)  estimate;  that  is,  the  minimiz- 
ing function  is  quadratically  normed  with  respect  to  the  inverse  of 
the  covariance  matrix  of  the  observation  errors . Precision  is  measured 
in  terms  of  a covariance  matrix. 


Definition.  Suppose  that  a minimum  variance  estimate  of  a set 


of  q parameters  X^  - jx-]_,  Xp,  . ..,  xc  ) has  been  obtained.  Then  we 


"til 

shall  call  the  i parameter  of  X an  uncorrelated  estimate  iff  x^  is 
uncorrelated  with  each  of  the  other  q - 1 parameters. 


Theorem  of  Uncorrelated  Estimation/  Let  X^  7 

be  a minimum  variance  estimate  for  a set  of  q parameters,  and  let 
y2>  •••  > yq j be  a q - dimensional  vector.  A parameter 
x^  is  an  uncorrelated  estimate  if  there  exists  a nonsingular  linear 
transformation  N such  that 
i ) Y - NX  , 


»Xj_,  x 2 } 


where 


N - 


— = 


I M 

o 1 

mi 

m2 


,1 


} 


m. 


q - 1 


s \ 
xi 


V 


xiX! 

3X2Xi 


xqxi 


and  ii ) at  least  one  of  the  transformed  parameters  in  Y is  uncor- 
related with  each  of  the  other  parameters  in  Y. 


Proof . The  proof  will  be  by  construction,  first  examining  the 
solution  for  two  parameters  and  then  for  q parameters.  There  will  be 
no  loss  in  generality  in  assuming  the  ith  uncorrelated  estimate  param- 
eter  is  in  the  q position,  and  will  be  designated  x . The  minimum 
variance  estimate  will  also  be  called  the  best  estimate.  First,  suppose 


1+8 


we  have  obtained  the  best  estimate  for  two  parameters,  and  Xg,  and 
their  associated  covariance  matrix 


We  want  to  find  the  condition  utilizing  the  covariance  terms  such  that 
the  resulting  parameters  are  uncorrelated,  and  such  that  at  least  one 
parameter  does  not  undergo  transformation.  Thus,  make  the  following 
nonsingular  linear  transformation  to  the  new  parameters  y^,  y2; 


choose  m such  that  y^  is  uncorrelated  with  yg  = x2-  That  is,  the  co- 
variance  of  y^  and  yg  is  zero: 


(10) 


or 


Y - LX  ; 


s - 0 


where 


8 


- E(Xi  4-  mx2  - Ex1  - mEx2)(x2  - Ex^ 

- E - Ex1)  4-  m(x2  - Ex2)  (xg  - Ex2) 

= E(X;l  - Ex1)(x2  - Exg)  + mE(x2  - Ex2)2 


Thus,  setting  s — 0 we  have 
*1*2 


4-  m s 


2 


s 


X. 


2 


49 


s. 


or 


m — • 


xlx2 

y. 


Since  the  values  of  the  covariance  elements  are  known,  the  value  of  m 
is  determined,  and  substituted  into  (10).  We  solve  (10)  for  the  orig- 
inal parameters  x^,  Xg* 

X - L_1  Y , 


“ 

-i 

- 

-1 

1 

m 

i 

— m 

where  L r 

_0 

ij 

— 

_0 

1_ 

and  substitute  L ^ Y for  X in  the  original  multivariate  regression 
analysis.  Our  minimizing  function  is  now  to  be  minimized  with  re- 
spect to  the  new  variables  y^,  y , which  are  now  uncorrelated.  Since 
at  least  one  of  the  new  parameters,  in  this  case  y^,  was  not  trans- 
formed, upon  obtaining  the  best  estimate  for  y^  and  y^,  the  value  of 
yg  - is  an  uncorrelated  estimate  and  will  be  fixed.  The  minimi- 
zation problem  has  now  been  reduced  one  dimension,  and  the  remain- 
ing parameter  X]_  is  solved  for  in  a new  regression  analysis,  yield- 

( ' ) (') 

mg  an  uncorrelated  estimate,  x^  ' , x£  ' for  the  parameters.  The 
whole  process  may  now  be  repeated,  where  the  new  value  of  m would  be 
obtained  from  the  new  elements  of  the  covariance  matrix  of  x|  ^ , x ^ ^ • 

8 , (’)„(' ) 


m 


(')  - 


X1  x2 
2 

SxO 

2 


designating  the  new  parameters  by  y|  \ yj^'). 
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We  shall  now  generalize  the  above  procedure  to  include  more 
than  two  parameters . 

Suppose  we  have  solved  a multivariate  regression  problem  for 
q parameters  x^,  x^,  . ..,  x^  together  with  the  associated  covariance 
matrix 


2 

Sxl 

Sx  V 
lx2 

Sxlx2  ' ‘ * 
s2 

8Xg  • * * 

1 

l s 

1 xlxq 

i 

i Sx2xq 

II 

* 

cal 

• 

• 

• 

1 

1 

• 

Sxlxq 

6x-x  • • • 

x2xq 

1 s2 

: x*  J 

the  matrix  as  shown 

: 

[ 

S - 
-x 

-11 

-^21 

-12 

S ' 

-22- 

_ 2 


where  has  dimension  (q-l)x(q-l)  and  - s 


-22  "q 

Make  a nonsingular  linear  transformation  to  a new  set  of  param- 


eters ylf  yg,  ...  , y : 


H* 

II 

I M 

k| 

y 

0 1 

X 

L qJ 

L qJ 

(11) 


or 


Y - NX 


the 


where  Y^  - j y^  yg,  . yq_]_|  , X^  - <jx1,  xg,  ...  , x j>  , I is 
identity  matrix  of  dimensions  (q  - l)  x (q  - l),  M^  - jm^,  mg,...,  mq 

and  0 has  dimensions  lx(q  - l).  We  want  to  choose  M so  that  the  com- 
ponents of  Y^  are  each  uncorrelated  with  yq.  Thus  the  (q  - l)xl 


51 


vector  M must  satisfy 


E(Y  - EY,)(y  - Ey  ) = 0 

-L  j-  <1  q 

Using  equation  (11 ) gives 

0 - E(XX  + Mxq  - - MExq)(xq  - Ex^) 

- e[  (X.  - EX.  ) + M(x  - Ex  )'  (x  - Ex  ) 

L x ~ q qjq  q 

= 3(3^  - H^Xx^  - Btq)  + ME(xq  - Exq)2 

= 512+«^2  • 


Then, 


M 


_ -12 
■‘T 

q 


since  S — s , a scalar.  Expanding,  gives 


mi 


( 


m 

r2  =- 


m 


q-i 


3jc2xq 


"Xq-lXq 


\ 


/ 


To  apply  this,  first  solve  equation  (ll)  for  X: 

X = N_1Y  , 


where 


N"1  — 


I -M 

0 1 


T-l, 


Substitute  N Y for  X in  the  original  regression  problem,  and  minimize 
with  respect  to  the  new  variables  YT=|y1,  yg,  ...  , y^  j>  , where  yx, 
Vo>  •••  > yn  i are  each  uncorrelated  with  y . Since  also  y was  not 

c.  q-x  q * q 

transformed,  the  value  of  y^  r:  x^  is  an  uncorrelated  estimate  when  a 
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best  estimate  is  obtained  for  Y,  the  vector  of  parameters.  This  value 
of  Xq  is  then  fixed,  which  reduces  the  minimization  of  the  remaining 
X^  - vector  elements  to  a (q  - l)  - dimensional  problem.  Designate 


this  new  solution  by 


X 1 -{xl'  x2'  *'*  > xq-l'  xq } 


where  only  x^  is  an  uncorrelated  estimate.  The  associated  (q  - l) 
x(q  - l)  - dimensional  covariance  matrix  is  written 


S = 
“*1 


D I 

xi 


S I , 

xlx2 


*x2 


s i , 
xlxq-l 


* ‘ I sx£x4_i 


-I 


xlxq-l 


s i . 

*2*4-1 


. . . i s i 

> xq-i 


where  the  partitioned  submatrices  are  given  by  the  array 


S'  3 
~x. 


-11 

-12 

s' 

s' 

“21 

“22 

where  s'  has  dimension  (q  - 2) x (q  - 2)  and  s'  = s^ 

-11  -22  x' 

q-1 


set  of  (q  - l)  parameters  y.J , y^,  ...  , y^_1  : 


As  before,  make  a nonsingular  linear  transformation  to  a new 

1'  J2 

I M Xg 

0 1 


l2 

Vl 


q-1 


(12) 


or 


= N'X^ 
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The  vector  M must  be  chosen  so  that  the  components  of  Yg  are  each  un- 
correlated with  yq_]_.  Following  an  identical  procedure  as  for  M,  we 
obtain  the  following  expression  for  M* : 


M - 


S 

-12 


<1-1 


2 

, where  S1  ~ s , 

-22 


q-1 


or 


/ • > 

ml 

“2 

. 

ii 

i 

ro 

( s ^ 

sx'x' 

xlxq-l 

sx'x' 

2 ^ 

sx’  n 

• 

^ mq-2, 

q-i 

V J 

Solve  for  N*  ^ from  (12 ) and  substitute  n'  _1  Y^  for  in 

the  (q  - l)  - dimensional  regression  problem.  Minimize  with  respect 

to  the  variables  Y^y'  y’,  ...  , y’  > , where  y*  ...,  y»  are 

1 1 2 q-lJ  1 q-2 

uncorrelated  with  y’  ..  . The  value  of  y1  =;  x’  is  an  uncorrelated 

<1-1  q-1  q_i 

estimate  when  a best  estimate  is  obtained  for  y\  Fix  this  value  of 
Then  minimize  and  obtain  a solution  for  the  remaining  - ele- 
ments. The  new  complete  solution  is  now  designated 

r*  f T 


— [ y’  ' x1  1 

r i'  x 2* 


x"  . X*  , X 
q-2  q-1  q 


where  now  x and  x'  are  uncorrelated  estimates. 

<1  q-1 

This  process  is  repeated  until  all  q parameters  are  uncorre- 
lated estimates,  resulting  in  a final  solution  vector. 


X(9-Dj  T = j (q-1)  (q-2) 

) (1  * x2  > •••> 


Xq-2»  xq-l'  xq 


5^ 


This  process  requires  a total  of  2q  - 1 regression  solutions  and  q - 1 
calculations  of  covariance  matrices.  The  entire  process  may  now  be 
repeated,  particularly  if  at  any  point  the  problem  vas  linearized  by 
the  usual  Taylor  series  expansion,  until  changes  in  the  final  solution 
vector  are  minimal  according  to  some  predetermined  criteria.  This 
completes  the  constructive  proof.  f 


We  shall  suppose  that  the  space  S governing  a set  of  vectors 
is  a complete  quadratic  finite  dimensional  pre-Hilbert  space.  Suppose 
that  ve  have  an  (mxn)  matrix  A of  rank  r=>0.  We  can  then  factor  A 
into  two  matrices  of  rank  r as  follows: 


The  columns  of  B form  a basis  for  the  column  space  generated  by  the 
column  vectors  of  A,  and  C is  chosen  to  satisfy  (13).  Each  mxl 
column  vector  of  B is  quadrat ically  normed  with  respect  to  an  mxm 
symmetric  positive  definite  matrix  M and  each  nxl  row  vector  of  C is 
quadratically  normed  with  respect  to  an  nxn  symmetric  positive  defi- 
nite matrix  N.  Then  we  define  the  quadrat 1 c pseudo inverse  of  A with 
respect  to  M and  N as 


Quadratic  Pseudoinverse 


A - B C . 
mxn  mxr  rxn 


(13) 


vm,n; 

If  A = 0,  then  we  define  0*~  0^ . 


(14) 


The  following  theorems  all  refer  to  matrices  and  vectors  as 


55 


defined  above . 


Theorem  3 • The  quadratic  pseudoinverse  of  an  mxn  matrix  A 
(with  respect  to  M and  N)  satisfies  the  property  of  a generalized 


inverse 


8 


AA , . A 3 A 

— (M,N)  ~ “ 


(15) 


Proof.  AA,„  XA  ~ BC 

— (m,n)-  — 


»“1CT(OT“1CT)'1(BT?ffi)‘1BrM 


BC 


= B(OT-V)(CT-V)'hBim)''L(BTMB)C 

= BIIC 
= A . 


■1  T,-l.  T. 


Theorem  k.  If  A is  nonsingular,  then  a7  ,=  A~^. 
- 6 “(M,N)  - 


Proof . We  know  that 


A_1A=AA‘1;=I, 


(16) 


since  A is  nonsingular  and  has  an  inverse.  Thus,  applying  A"1  to 
both  sides  of  (15)  on  the  left  gives 

A-1(AA*  Afc A-1A  , 

_ ~“(m,n)  ~ 


(A  A) (A*  A)=A_iA  , 

(M,N) 

which  implies  A^  ^A-I^  , aright  identity  matrix. 

Now,  apply  A ^ to  both  sides  of  (15)  on  the  right,  giving 


(IT) 


(18) 
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which  implies  AA^  = IL  , & left  identity  matrix.  From  (l6),  (17) 

and  (18)  we  have  I L-I-I  I implies  I = I implies  A*  „SA  = AA* 

- ^ -X  HR  F “(M,N)-  — (M,N) 


= A"1. 


implies  A. 

~(m,n)  - 


I 


Theorem  5.  (A*  )*  —A  . 

~(M,N)  (M,N)  - 


Proof.  Write  (l4)  in  the  form 


where 


Then 


%,*r- 


Q = n~1ct(cn"1ct)'1  , P:r  (bt^)_1b?m 

QT—  (CH_1CT)_1OT”1  and  PT= 


By  definition  of  the  quadratic  pseudoinverse  of  A*  with  respect  to 
M and  N,  we  form 

Direct  substitution  of  the  above  equalities  for  Q and  P yields 

(-(M,N)}(M,N)r  — • | 


Theorem  6.  Let  the  rank  of  the  mxn  matrix  A be  n.  Then 
L(M) 


hju\  -(atma)_1atm  . 


Proof . If  r-  n,  then  from  equation  (13), 

A — B since  C — I . 
mxn  mxn  nxn  — 
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Direct  substitution  into  the  definition  (l4)  yields 
A*  = N“1I(IN"1I)"1(ATMA)“1ATM. 


Since  the  quadratic  pseudoinverse  is  independent  of  N when  the  rank 
r of  A is  n we  may  designate  the  dependence  on  M alone  by  A*^ y Thus 


A*  = (AtMA)_1ATM  . 
— CM; — 


I 


Theorem  7 • Let  the  rank  of  the  mxn  matrix  A be  m.  Then 
A*n  j - N-V  (AN  - 1AT)‘1. 


Proof . If  rrm,  then  from  equation  (l3)> 

A — C since  B =1  . 

mxn  mxn  mxm 

Substitute  into  equation  (l4)  to  give 

A*  = N'1CT(CN'1CT)_1(IMI)-1IM. 
~(M,N)  - - — - — 


Since  the  quadratic  pseudoinverse  is  independent  of  M when  the  rank  r 


of  A is  m we  may  designate  the  dependence  on  N alone  by  Then, 

a*n)=n'1at(m-1at)-1.  I 


Theorem  8.  The  quadratic  pseudoinverse  of  a product  is  the 
product  of  the  quadratic  pseudoinverses  in  reverse  order,  all  taken 
with  respect  to  M and  N. 


Proof . It  is  sufficient  to  consider  the  product  in  equation 
(13),  since  the  extension  is  immediate.  Thus 
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A,  . - (BC) , “ C B 

-(M,N)  V~  (M,N)  -(H)-  (M) 

by  applying  the  tvo  previous  theorems . 


Quadratic  Pseudoinverse  Estimation 
This  estimation  is  an  extension  of  that  developed  in  the  sec- 
tion on  Multivariate  Regression  Analysis.  In  that  section  the  nxq 
matrix  A had  rank  q,  q-=:n.  This  section  generalizes  the  nxq  matrix 
A by  assuming  its  rank  is  r5qcn(  i.e.,  only  r column  vectors  of  A 
are  linearly  independent . 

As  in  the  previous  sections  ve  shall  suppose  that  the  space  S 

governing  the  vectors  is  a complete  quadratic  finite  dimensional  pre- 

Hilbert  space  and  that  N and  Q are  symmetric  positive  definite 

nvn  qxq 

matrices . 


Theorem  9.  Suppose  that  the  nxq  matrix  A has  rank  r such 
that  r=  q<n.  Suppose  that  we  can  factor  A into  tvo  matrices,  each 
of  rank  r,  such  that 


A - B C 
nxq  n*r  rxq 


(19) 


where  the  columns  of  B form  a basis  for  the  column  space  generated 


by  the  column  vectors  of  A,  and  C is  chosen  to  satisfy  (19).  Sup- 
pose that  each  nxl  column  vector  of  B is  quadratically  normed  with 
respect  to  an  nxn  symmetric  positive  definite  matrix  N and  each 
qx  1 row  vector  of  C is  quadratically  normed  with  respect  to  a qxq 
symmetric  positive  definite  matrix  Q.  Then  the  parameter  estimate 
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A 

x for  the  absolute  minimum  of  the  quadratic  minimizing  function 

M(x)  =||  y - Ax  '|2  , 

where  y is  the  nxl  vector  of  measurements  and  x is  the  qxl  vector 
of  parameters,  is  given  by 

x = Q"1CT(CQ"1CT)"1(^^)’1BTNy  . 

Proof . Since  A — B C,  we  can  write 
M(x)  — 1 1 y - Bz  ||  2 

where  z — Cx  and  now  the  n x r matrix  B is  of  rank  r,  z is  an  rx  1 vec- 
tor, and  the  conditions  of  the  first  theorem  in  the  Multivariate  Re- 
gression Analysis  section  are  satisfied.  Thus  the  estimate  of  z is 

z = (BTNB)'1BrNy  . (20 ) 

Since  the  rxq  matrix  C is  composed  of  known  elements, 

z - Cx  . (21) 

Combining  (20 ) and  (21)  and  solving  for  x by  taking  the  quadratic  pseu- 

doinverse of  C with  respect  to  the  qxq  symmetric  positive  definite 
matrix  Q gives 

* = £*Q)(BTNB)-1BTNy  • 

Since  C is  an  rxq  matrix  of  rank  r, 

£(q)  “ Q“1ct(cq_1ct)'‘1  . 

Hence,  x = Q“1CT(^"1CT)'1(BTI^)‘1BTNy  . (22 )| 

By  definition,  we  see  that  the  coefficient  of  y in  equation 
(22)  is  the  quadratic  pseudoinverse  of  A with  respect  to  N and  Q. 


6o 


That  is 


Constrained  Estimation 


When  an  a priori  estimate  of  the  variance  of  a parameter  is 
known,  it  becomes  desirable  to  constrain  the  variance  by  this  a priori 
estimate.  Experience  has  indicated,  according  to  D.  Brown,  9 that  the 
inverse  of  the  variance -covariance  matrix,  given  by 


in  a minimum  variance  estimation,  turns  out  to  be  nearly  singular.  This 
is  attributed  to  linear  dependence  of  two  or  more  of  the  parameters, 
and  constraining  the  variances  causes  the  resulting  matrix  to  be  non- 
singular. The  following  definitions  and  theorems  examine  this  new  con- 
cept. We  assume  that  the  space  S governing  the  vectors  is  a complete 
quadratic  finite  dimensional  pre -Hilbert  space  and  that  the  weighting 
matrices  are  symmetric  positive  definite. 


P = ATS  ~1A 

•J 


Definition . Suppose  that  a minimum  variance  estimate  of  a set 


rn 

of  q parameters , x 1 
call  the  i^*1  parameter  of 


parameter  of  x a constrained  estimate  iff,  by  assuming  the 


obtained.  Then  we  shall 


covariances  of  the  parameters  equal  to  0,  the  variance  of  the  ith  param- 
eter is  bounded  above  by  the  a priori  variance . 


The  following  theorem  is  based  on  the  derivation  by  D.  Brown. 


minimum  variance  estimate  for  a set  of  q parameters.  A parameter  x± 
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is  a constrained  estimate  if  the  regression  analysis  yields  the  pa- 
rameter estimation 

x =(FT^)"1rTSz  , 
where  the  weighting  matrix  S - 


-T  J- 

X - < X.  , X, 


R 

0 

r n 

A 

- ~] 
y 

0 

Q 

-a  — —y  > z= 

I 

, z~ 

aJ 

and 


1*  2*  •••>  , 


i/sf 


Q = 


1 


i hi 


l/s£ 

X 


2 


il 


such  that  s~  is  the  a priori  variance  of  the  ith  parameter  and  x 
xi  i 

is  the  a priori  value  of  the  1th  parameter. 


The  matrices  A and  R and  the  measurement  vector  y are  the  same 
as  in  the  section  on  Multivariate  Regression  Analysis. 


Proof.  To  the  vector  y - Ax  in  the  quadratic  minimizing  func- 
tion add  the  non-zero  vector 

x - I x 

where  x is  the  vector  of  a priori  values  of  the  parameters: 

y - A x 


x - I x 


- z - F x 


(23) 


where  F r: 

(n+qjxq 


z - 
(n-i-q)xl 


y] 

*J 


The  covariance  matrix  of  the  combined  vector  is 
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S'1  = 


where  S — 


S 

nx& 

0 


s£ 

X, 


0 


sr 


Sx 

qxq 


s£ 

Xr 


(2k) 


, the  a priori  covariance 


matrix  of  the  parameters . 

From  equations  (23)  and  (24),  ve  can  form  the  quadratic  mini- 
mizing function 


S(*>=||z  - £*11  s > 


provided  that  S ^ has  an  inverse . The  minimum  variance  solution  is 

/n  rp  _ 1 rp 

x = (f SF)  F S z , 

where  the  covariance  matrix  of  the  parameters  is 


£U  = (FTSF)_1 
~x 


s'1 

-y 

o 


rp 

S.  = I A1SV  A + S 

r ~y  ” "x  j 


0 

§.  "j 


-1 


J 


, where  S- 


ATS-1A  J -1  . 

--y-J 


Let  us  designate  the  ith  diagonal  element  of  S_1  by  b^.  If  we  assume 


that  the  covariances  of  the  parameters  are  equal  to  zero,  then  the  i 

diagonal  element  of  SZ1  is  given  by  l/s£  . Then 

2 Xi 
1 b 1 

= i + 


th 


sx. 
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bo  that  s|  = s£  , 

i xi 

t h 

that  is,  the  variance  of  the  i parameter  is  bounded  above  by  the  a 
priori  variance . jjj 


Functionally  Dependent  Constrained  Estimation 
In  the  case  that  there  is  a knovn  functional  dependence  among 
the  errors  in  the  parameters,  ve  can  derive  a nev  a priori  covariance 
matrix . Suppose  ve  can  vrite  the  vector  difference  as 

dx  - X—  x , 

vhere  x is  the  vector  of  a priori  values  of  the  parameters  and  x is  the 
true  vector  of  parameters.  But  if  a functional  relation  exists  among 
these  errors,  then  there  exists  a singular  matrix  G such  that 

Gdx  - x - x . 

This  leads  to  the  following 


Theorem  10.  Suppose  a knovn  functional  relation  exists  among 
the  errors  in  the  q parameters  of  a minimum  variance  estimate  and 
that  there  exists  a singular  matrix  G defined  by 

Gdx  - x - x , 

vhere  x is  the  a priori  parameter  vector  and  x is  the  vector  of  true 
values.  Then  the  a priori  covariance  matrix  of  the  parameters  is 
given  by 


S_ 


G S G 
dx- 


> 


vhere  S is  a diagonal  matrix  of  the  a priori  variances  of  the  pa- 


rameters . 
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Proof.  This  is  immediate  from  the  section  on  derivation  of 
covariance . | 

We  may  now  combine  the  above  theorem  with  the  Theorem  of  Con- 
strained Estimation  to  produce  the  following 


Theorem  of  Functionally  Dependent  Constrained  Estimation. 
Suppose  that  a known  functional  relation  exists  among  the  errors  in 

the  q parameters  of  a minimum  variance  estimate,  given  by 

^T  _/  l 

x — ) x^  t Xg , . . . , x > 

such  that  there  exists  a singular  matrix  G defined  by 
Gdx  — x - x , 

where  x is  the  a priori  parameter  vector,  x is  the  vector  of  true 
values  and  dx  is  the  vector  of  errors.  Let  us  call  the  ith  param- 
eter of  x a functionally  dependent  constrained  estimate  iff,  by 

assuming  the  covariances  of  the  parameters  equal  to  zero,  the  vari- 
"t»h 

ance  of  the  i parameter  is  bounded  above  by  a function  of  the  ele- 
ments from  the  a priori  covariance  matrix  of  the  parameters ? Then 
the  parameter  x^  is  a functionally  dependent  constrained  estimate  if 
the  minimum  variance  estimate  is  given  by 

x =(FTSF)"1FTS  z 


where 


The  weighting  matrix  S is  given  by 


L°  aj 


* 


As  given  in  Theorem  10. 
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where  R 1=  £L,  and  Q"1—  S_  = G S,  GT . 

— -v  ~5c dx— 

The  matrices  A and  R and  the  measurement  vector  y are  the  same  as 
in  the  Theorem  of  Constrained  Estimation.  S-  is  the  a priori  co- 
variance  matrix  of  the  parameters  and  is  a diagonal  matrix  of 
the  a priori  variances  of  the  parameters . 


Proof . The  proof  is  identical  to  that  of  the  Theorem  of  Con- 


strained Estimation  except  that  there  must  be  replaced  by 

x 

S = Q S GT  . 
dx 

Then  the  covariance  matrix  of  the  parameters  is  written 


S.i'aYaHGS  gt)_1 

-x y - - -dx~ 


-1 


The  variance,  due  to  the  measurements  only,  of  the  i^*1  parameter  is 
given  by  the  ith  diagonal  term  of 


S = 


-1 


rth 


Designate  the  ix  diagonal  element  of  S by  b^.  If  we  assume  the  co- 
variances  of  the  parameters  are  equal  to  zero,  then  the  ith  diagonal 
element  of  is  given  by  l/s?.  Designate  the  i^*1  diagonal  element 


—x 


of  ^ ) 1 > where  the  ith  diagonal  element  of  G S.  GT  is 


■f  "h 

the  a priori  variance  of  the  i parameter.  Then, 


S2 

x 


1 = bI 


>i  > 


which  implies  s^  ^ l/g^  . 


I 


66 


Derivation  of  Covariance 

According  to  definitions  in  Wilks,  reference  (}D),  the  follow- 
ing statistical  properties  of  random  variables  may  be  obtained.  The 
expected,  or  mean,  value  of  a random  variable  u defined  by  equation 
(a)  page  29,  of  ref . (ID)  will  be  designated  by  the  symbol  E(u).  The  co- 
variance  of  two  random  variables  u,  v is  given  by 


s = E 
uv 


u - E(u) 


E(v  ) 


= E juv  - vE(u)  - uE(v)  + E(u)E(v)\ 

- E (uv)  - E(u)E(v)  . 

If  u and  V are  independent  of  each  other,  then  the  covariance  is  zero. 
If  u — v,  the  covariance  reduces  to  variance 

rE(u2)  - !E(u)]2. 


Theorem  12.  Suppose  vj  ,(j  — 1,  ...,  n)  are  n random  variables 
such  that 

n 

} • • • y HI ) > 

j=i  J J 

where  the  a^j  are  constants . In  matrix  form  this  may  be  written 


/ \ 

; ui 

all 

a12 

• ' ’ aln 

fvi' 

u2 

a21 

a22 

...  a2n 

- v2 

, 

{ . 

\= 

• 

/ 

\ • 

• 

• 

• 

• 

• 

1 • 

! U 

l V 

a -i 

ml 

am2 

• • • awy, 
mn 

or  U = A V 
m*l  mxn  nxl 


Define  the  following  covariance  matrices : 
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s = 

-u 
m xm 


u. 


S 8 

^Ug  Ug 


Vm 


UlUg 

2 


u2um 


b . . . s 

u0u  u 

2 m m 


S : 
—v 
n Xn 


Then 


V-i  v, 


lv2 


v v 

1 n 


Vlv2 


Vlvn 


• • • s 


^2vn 


Vn 


. 8 


s = A S A . 

-u v — 


Proof • Fr°®  the  linear  properties  of  the  expected  value  oper- 


ator. 


E(\  u-^-E 


(£ 

,1  = 1 
n n 


aki  vi 


)(  ^ aij  Vj) 

J=1  J 


=*(£  L 


J=1  i— 1 


aki  al j vi  vj ) 


E(u  ) E(u  ) 
k 1 


J = 1 


- \i  aij  E(vl  v,  > 

i=l  J 


= E(£ 


i = l 


a,  v 
ki  i 


) E(  Y 

J = 1 


1J 


v 


~n 

n 

E 

£ alj  E(vi} 

i = 1 

J=1 
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n n 

= L L aki  aii  E(vi^  E(y<) 
j-i  1=1 


Vl  ^E(uk  V " E(V  E(ul} 


n n 


“ ^ £.  V \j  E(vi  V ' ^ ? °M  alJE(vl)E(vj) 

- x 1 ~ x j = 1 i - 1 


j 


n n 


= £ £ \i  au 

J = 1 i — 1 


E(vi  vj^  “ E^vi^  EM 


n n 


\\=Z  £ 


a . a.  , s 
ki  lj  v.v 

J = 1 i-1 


Direct  calculation  of  the  kl -element  of  A A1  yields  the  same  result 

as  above  for  s . Therefore 
ukul 


S “AS  A 
"hi v “ 


as  desired. 


I 


12 

Definition.  Two  random  variables  u,  v are  uncorrelated  iff 
E (uv)  = E(u)  E(v)  , u±v 

If,  in  addition 

E (uv)  — 0 

then  the  two  random  variables  are  orthogonal . From  this  we  see  that 
orthogonal  random  variables  with  zero  means  are  uncorrelated. 


12  „ 


Definition.  Two  random  variables  are  independent  iff  they  are 
orthogonal  random  variables  with  zero  means  (uncorrelated)  and  have 
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normal  (Gaussian)  distributions. 

Corollary  2 . Suppose  the  v , (1  = 1,  . ..,  n),  are  n independent 
random  variables.  That  is,  they  are  normally  distributed  and 
E(v!  vj)=E(vi)E(vJ)=0,  i^J,  (i,  J=l,  ...,  n). 

Then 

S -AS  AT 
-u  v “ 

where  S r 
V 


Probability  Ellipsoid  of  Error 3 A3 
The  elements  of  the  variance -covariance  matrix  may  be  used  to 
describe  an  error  ellipsoid,  which  is  centered  about  the  best  estimate 
of  position  (x,  y1,  £),  and  within  which  there  is  a predetermined  proba- 
bility of  the  true  position. 

Let  the  variance -covariance  matrix  S of  x,  y,  z be  given  by 


4 

sx2 

X 

xy 

s~~ 

8? 

xy 

y 

yz 

S^/s 

xz 

A 

yz 

t 

where  S/^  is  an  estimate  of  the  covariance  between  the  random  variables 
u and  v. 


TO 


If  we  assume  that  the  distribution  is  normal  then  the  equation 
of  the  ellipsoid  is  given  by 


G(x,y,z)  = 


x - x,  y - y,  z-z 


’ll 


12 


s 


12 

322 


13  23 


“l3 

323 

533 


A 

X - X 


y-y 

A 

z-z 

L j 


- cr 


(25) 


where  b±.  is  the  cofactor  of  the  ith  row,  jth  column  element  in  S, 
divided  by  the  determinant  value  of  S,  (x,  y,  z)  is  the  best  estimate 
of  position,  and  C is  a constant. 

The  probability  p that  the  best  estimate  of  position  (x,  y,  z) 
will  fall  within  the  ellipsoid  corresponding  to  a particular  value  of 


C is 


.2  , -*-2 
t e 


■td/2 


dt 


The  following  table  shows  some  approximate  values  of  C corresponding  to 


P percent  probability  ellipsoids  (P~100  p): 


p 

25 

50 

75 

82.8 

90 

95 

99 

c 

1.101 

1.538 

2.027 

2.236 

2.500 

2.795 

3-368 

Without  loss  of  generality  we  may  set  (x,  y,  z)  to  zero.  Then 
the  equation  of  the  ellipsoid  is  centered  at  (0,  0,  0),  and  may  be  re- 
duced to  diagonal  form  by  an  appropriate  orthogonal  transformation, 
from  the  (x,  y,  z)  coordinate  system  to  a new  system  (x,  y,  z),  where 
the  axes  of  the  ellipsoid  lie  along  the  axes  of  the  new  system.  In 
this  new  system  the  equation  of  the  ellipsoid  will  appear  as 
a^  + agy2  + a^z2  - C2  , 

or,  in  matrix  notation  corresponding  to  equation  (25)  it  will  be 


(26) 
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x,  y,  z 


a]_  0 0 

0 

0 0 


a2  ° 


x 

y 

z 


= C 


The  elements  a^,  a a^  of  the  diagonal  matrix  above  are  the  eigen- 


values of 


LsiJJ 


, and  are  obtained  from  the  eigenequation 


al 


v o 


where  1^  is  the  identity  matrix. 

We  may  reduce  equation  (26)  to  normal  form  by  dividing  through 


by  C and  rearranging: 


+ 


— 2 — o 

JL 4. 5_ 


1 . 


(C2/a  ) (C2/a  ) (C2/a  ) 


i'  vw  ' “2y  ’~y 

Thus,  the  lengths  of  the  semi  - axes  of  the  ellipsoid  of  error  in  the 
transformed  coordinate  system  (x,  y,  z)  are  VC/V  , (l-l,  2,  3). 

The  eigenvectors  of  the  semi  - axes  of  the  ellipsoid  may  be  obtained 
in  terms  of  the  direction  cosines,  with  respect  to  the  original  (x,y,z) 
coordinate  system,  by  solving  the  sets  of  simultaneous  equations, 

(1=1,  2,  3)  , 


jk 


T 

where  - 

to  each  a^ . 


V.q,  V^2 , Vi^  , for  the  eigenvectors  V corresponding 


x,  y 


- ai 

Vi  ' 

vi3^ 

, for 

ellipsoid 

in 

“811 

S12 

S12 

8 22 

- C 


(27) 
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vhich  is  the  equation  of  an  ellipse  in  general  form,  centered  at  the 
best  estimate  of  position  (0,  0).  The  eigenvalues  and  eigenvectors 
for  the  ellipse  are  obtained  by  methods  analogous  to  those  for  the 
ellipsoid.  If  b^  and  bg  are  the  eigenvalues  of  the  square  matrix  in 
equation  (27),  then  the  lengths  of  the  semi  - axes  of  the  ellipse  are 

V^/b.  , (1  = 1,  2). 


If  we  assume  that  the  distribution  is  normal,  then  the  proba- 
bility p corresponding  to  a value  of  C is  given  by 

- , -C2/2 

p — 1 - e ' 

The  following  table  lists  approximate  values  of  C for  probability 
ellipses  (P~100  p). 


p $ 

25 

50 

75 

86.466 

90 

95 

99 

c 

• 7585 

1.177 

1 .665 

2.000 

2.146 

2.448 

3.035 

Matrix  Calculus 

An  n-dimensional  vector  x is  written  as  a one-columned  matrix: 

/ „ \ 


x = 


a2 


x 

\ n/ 

Whenever  possible,  the  elements  of  a vector  will  be  contained  within 
braces  as  above.  A rectangular  mxn  matrix  of  m rows  and  n columns 
will  be  underlined,  e.g.  A.  If  necessary,  the  order  will  be  designated 
beneath  the  matrix.  The  elements  of  a rectangular  matrix  will  be 
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contained  within  square  brackets.  Thus,  if  a^  are  the  elements  of  a 


matrix  A of  order  mxn,  we  shall  write 


A - 
mxn 


Certain  matrix  operations  are  denoted  by  superscripts.  Thus 
T 

A , transpose  of  A 
nx  m mxn 

A ^ , inverse  of  A , defined  only  for  square  ma- 
mx m m x m 

trices . 

Define  the  matrix  partial  derivative  operator  with  respect  to 
the  q variables  q^,  q^,  ...,  qn,  as  a vector 

1 b ' 


n Xl 


< : 

• 

b 

ba 

\ / 

Applied  to  the  transpose  of  an  mxl  vector  x whose  elements  are  func- 
tions of  q we  have 

( b \ 


Vfl  T - 
q x ~ 

nx 1 1 Xm 


) { xl'  * ‘ ' *01  } 


J>_ 

V \ J 


7k 


dq 

T1 

If  x is  a scalar  the  above  reduces  to 


6: 

S] 


m 


)(1n 


3 B 
n x m 


V, 

n x 1 


X 3 


V. 


d: 


bx 


b 
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and  if  V consists  of  a single  derivative,  say  — t we  have 
^ at  7 


V t 

lxl  Ixm 


,T  x <L 
dt 


m 


l 


dx-j 

dt 


dt 


^ T 

- y 


Inverse  of  Partitioned  Matrix 

Theorem  13  ■ Let  V be  an  nxn  matrix  and  partition  the  matrix 
as  follows : 


V = 


V-1  cc 


hi 

r x r 


I 


s x r 


*12 

r x s 


~22 
s x s 


, where  r+s-n. 


hi  - hshz  hi ] -1 


- hi  h.2 


-1, 


r 


-4g 


r-l 


2 hl[hl  - —12  —22  hi 


Y22  ■ ¥21  V1L  V^l 
^22  “ -21  -11  —12] 


Then 
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Proof.  Find  V-1  such  that  V V"1^!  . That  is 


Xn  Vd2 


L,  v 

-21  -22 


H12 


u u 

"21  “22 


n x n 

I 0 

rxr 

0 I 

S X s 


This  results  in  the  following  four  equations: 


a)  Xufiu+ViaS^i 
(1')  r*  r and  (2') 

b * -21  -11+  hs  -21  = 0 


-11 

—12 

where  V-  — 

u 

U 

-21 

-22 

hi  —12  "*"—12  b22 

= 0 

b>  hi  “12+^2  5*22 


: I 

S X s 


Solve  the  set  of  equations  (l1)  simultaneously: 


Hgi  ~ ” ^2  —21  —11  Prom  O an<l  substitute  into  a),  giving 


-11  -11  ‘ —12  -22  hi  hi  - i lmplles  hl= 

rxr  xx 


V - V V”'*’  V 
-11  -12  -22  -21 


Then  U — 
“21 


V-1  V 
“22  -21 


V - V V-1  V 
-11  -12  -22  -21 


-1 


In  a similar  fashion,  solving  equations  (2')  gives 

-1 


^22  " 


—12  ~ 


V - V V -1  V 
-22  ^21  -11  -12 


and 


v”1  V 

-11  -12 


v - v V-1  V _1 

^22  -21  -11  —12 


Corollary  3 • Let  and  be  null  matrices  in  the  above 


partitioning,  so  that. 


V = 


hi 

rx  r 
0 


-21 


^22 
S X s 


Then, 


-1 

1 

1 

V, 

0 

1 

-11 

0 

V'i 

“22 

1 

Cor°lla.ry  k.  Let  V be  a diagonal  matrix  with  elements  a 
the  real  scalar  field: 


V - 


Then, 


l/a-L  0 

1/ag 

l/a 
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CHAPTER  IV 


OPTIMIZATION  OF  SATELLITE  DOPPLER  NAVIGATION 


Introduction 


We  shall  apply  optimization  theory  to  the  problem  of  ship  nav- 
igation by  the  observation  of  radio  Doppler  signals  from  near  earth 
satellites.  Soon  after  the  first  man-made  satellite  was  placed  into" 
orbit,  Guier  and  Weiffenbach  of  The  Applied  Physics  Laboratory,  The 
Johns  Hopkins  University  (APL/jHU)  shoved  in  reference  (l4)  that  it  is 
possible  to  determine  the  orbit  of  a satellite  that  transmits  a fre- 
quency-stable continuous  wave  (CW)  radio  signal  by  measuring  the  Doppler 
shift  of  the  signal  on  ground  based  receivers.  In  1962  (ref.15) 
Kershner  and  Newton,  also  of  APL/JHU,  demonstrated  that,  having  deter- 
mined the  orbit,  the  technique  may  be  inverted  and  the  position  of  a 
receiver  on  the  earth  may  be  obtained  by  observing  the  Doppler  signal 
from  the  satellite.  Each  estimation  method  will  be  applied  to  the 
more  difficult  problem  of  solving  not  only  for  the  position,  but  also 
for  the  velocity  of  a receiver  on  a moving  ship  from  a single  15- min- 
ute pass  of  a satellite. 

The  basic  measurement  is  the  observed  frequency  of  the  CW 
transmitter  of  the  satellite . If  the  satellite  is  radiating  at  the 
frequency  f,  then  the  measured  Doppler  signal,  A f,  is  given  by 


f ] dA 
c i dt 


(1) 


77 


78 


where  A is  the  optical  path  length  between  satellite  and  receiver; 
is  the  time  rate  of  change  of  A ; 
f is  the  transmitted  frequency;  and 
c is  the  velocity  of  light. 

If  there  were  no  refraction  along  the  optical  path,  dA  /dt  would  be 
just  the  time  rate  of  change  of  range.  In  this  derivation  we  shall 
assume  no  refraction,  or  that  a refraction  correction  has  been  made. 
Designating  the  range  between  satellite  and  receiver  by  p,  we  may  re- 
place dA  /dt  in  equation  (l)  by  £. 

In  practice,  the  actual  Doppler  signal  observations  are  taken 
by  counting  a fixed  number  of  cycles  of  the  transmitted  signal,  measur- 
ing the  span  of  time  within  which  they  occur,  and  initiating  this  count 
at  intervals  of  1,  2,  or  k seconds  throughout  the  time  span  during 
which  usable  data  may  be  acquired.  A single  such  observation  of  the 
satellite  is  known  as  a "pass.  " 

Multivariate  Regression  Analysis  of  Satellite  Doppler  Navigation 
The  basic  technique  of  ship  navigation  by  the  Doppler  obser- 
vation of  satellites  is  to  obtain  the  best  estimate  of  the  ship's 
position  and  velocity,  which  results  in  a minimization  of  the  sums  of 
squares  of  the  difference  between  the  observed  Doppler  signal  and  the 
Doppler  signal  computed  from  the  positions  and  velocities  of  the  ship 
and  the  satellite.  This  statistic  also  produces  a probability  ellipse 
of  the  error  in  navigation.  The  difference  between  the  observed  and 
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the  computed  Doppler  is  the  residual  Doppler  signal. 

To  apply  the  theorem  of  multivariate  regression  analysis  to 
navigation  by  Doppler  observations  ve  shall  investigate  in  more  detail 
the  following:  l)  derivation  of  range  rate,  shoving  the  dependence  of 

the  Doppler  signal  on  the  positions  and  velocities  of  the  ship  and  the 
satellite;  2)  establishment  of  the  minimizing  function;  3)  lineariza- 
tion of  the  expression  for  range  rate,  making  the  problem  tractable; 
k)  minimization  of  the  minimizing  function,  resulting  in  the  set  of 
normal  equations  for  the  best  estimate  of  the  ship's  position  and  ve- 
locity. 

It  was  observed  above  that  the  Doppler  signal,  A f , is  directly 
proportional  to  the  range  rate,  p,  between  the  satellite  and  the  ship: 


If  we  let  capital  letters  indicate  the  actual  position  of  the  satellite 
and  lower  case  letters  the  actual  position  of  the  ship,  then  the  range 
p between  the  satellite  and  the  ship  is  given,  in  a rectangular  Carte- 
sian coordinate  system  by 


Differentiating  the  above  with  respect  to  time  and  dividing  both  sides 
by  p yields 


where  the  dot  indicates  time  differentiation.  The  above  equation  is 
valid  for  a receiver  on  a moving  ship. 

We  now  make  the  restriction  that  we  know  the  ship's  equation 


(!') 


p2  = (X  -x)2  + (Y  - y)2  + (Z  - z)2. 


P = i (x  " X)(x  - x)  + (Y  - y)(Y  - y)  + (Z  - z)(Z  - z)  , (2) 


8o 


of  motion.  For  our  purposes  we  can  assume  that  the  ship  moves  at  con- 
stant velocity.  Then  the  position  at  time  t may  be  found  in  terms  of 
the  position  (xq,  yQ,  zQ)  at  time  tQ  and  the  velocity. 


x 


'yW  y0Wsy  S (t  - tQ) 


(3) 


o 

\ j 


1 


Substitute  equation  (3)  into  (2)  and  then  (2)  into  (l1)  to  give  the 
following  for  each  i"*'*1  position  and  velocity  point  of  the  satellite 
at  time  t . : 


(Xi-x0-x(tj-t0))(X1-x)-KYi-yo-y(ti-to))(Yi-yH-(Zi-20-z(t1-t0))(z;1-z) 
[(Xi-xo-x(ti-to))2-t-  (Yi-y0-y(ti-t0))2+  (Zi-z0-z(ti-t0))2j  £ 

This  completes  task  l)  outlined  above.  (4) 

To  establish  the  minimizing  function,  first  we  must  state  the 
residual,  that  is,  the  difference  between  the  observed  and  computed 
Doppler  signals.  Let  the  observed  signal  vector  be  indicated  by  Af, 
and  the  computed  Doppler  signal  vector  by  Af’;  the  difference  is  then 
the  residual:  ( Af  - Af1).  Although  the  observed  Doppler  signal  is  re- 
ceived as  an  analog  signal,  that  is,  a continous  wave,  it  is  converted 
to  digital  (numerical)  form  at  discrete  time  intervals  t^,  and  the 
computed  Doppler  signal  is  calculated  from  equations  (4)  above  at  the 
same  times  t, . Thus,  for  a single  satellite  pass  there  are  n discrete 
residuals,  the  ith  residual  being  indicated  by  the  subscript  i.  Then 
applying  multivariate  regression  analysis,  we  assume  that  the  best 
estimate  of  the  ship's  position  and  velocity  (x^,  y^,  Zq,  x,  y,  z)mx'^ 
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can  be  found  by  minimizing  the  quadratic  minimizing  function 

S(x’)  =||  Af  - Af 'I]  ^ (5) 

where  R is  the  known  weighting  matrix,  Af  is  the  observation  vector 
and  Af'  is  the  vector  whose  elements  are  computed  from  equation  (4). 

In  general,  the  minimization  is  done  by  finding  the  partial  deriva- 
tives of  S with  respect  to  x£,  y£,  z£,  x}  yj  zj  setting  the  results 
equal  to  zero,  and  solving  for  (x^,  y£,  z'Qt  x\  yj  z ').  Due  to  the  non- 
linear character  of  equation  (4),  the  method  of  solution  for(x',  y'. 

x o 7 w o 7 

zo>  x}  y>  zj ) presents  problems  which  are  overcome  by  linearizing  the 
equation  so  that  the  theorem  on  multivariate  regression  analysis  can 
be  applied. 

The  linearization  is  achieved  by  expressing  the  i^*1  vector  ele- 
ment, A f£, . in  a Taylor  series  expansion  about  some  preliminary  (zeroth) 

estimate  of  the  ship's  position  and  velocity,  which  we  shall  indicate 

T 

by  replacing  the  prime  by  zero,  i.e.,  x°  = (x°,  y°,  z°,  x°,  y°,  i°). 

Then  Af^  is  the  theoretically  exact  Doppler  signal  computed  from 
equation  (4)  if  the  ship  is  actually  at  (x°,  y°,  z°,  k° , y°,  z°)  when 
the  satellite  is  at  its  ith  position.  The  expansion  of  Af^  about  the 
zeroth  estimate  shows  only  the  constant  and  linear  terms: 


Af 


i 


. o 

= Af±  + 

+ 


I 


0 6Af' 

^Af! 

Xq)  + by^y° 

- y0) + > , (zo  - *§) 

- zo 

bhfi 

dflfj 

- x°)  + — — (y* 

- y°)  + — -(z ' - z°) 

by' 

C)z ' 

+ higher  order  terms, 


(6) 
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where  each  partial  is  evaluated  at  the  zeroth  estimate,  and  is  computed 
for  n positions  of  the  satellite.  This  is  substituted  into  equation 
(5)  which  now  becomes 


S — 1 1 d ( flf ) - Adx|j‘ 


where 


d(  flf)  = 
nxl 


= flf  - Af°  - 


flf3 

fl  tr 


A f 


- flfj 

- 1st  2 


Af° 

n 


yo 


dx  n x ' - x° 
6x1 


o \ 


y' 


J 


A is  the  matrix  of  partial  derivatives  evaluated  at  the  zeroth  esti- 
mate . 

Apply  the  multivariate  regression  analysis  theorem  to  S to 
yield  the  parameter  estimate 


dx  = (/  RA)_1ATR  d(  Af). 


(7) 


LetS^  represent  the  covariance  matrix  corresponding  to  measurement 
errors,  let  £h  represent  the  covariance  matrix  corresponding  to  param- 


eter estimate  errors.  Then  the  weighting  matrix  R is  given  by 

,-l 


R — S 


'Af 


83 


equation  (7)  becomes 


dx-x'  - x°=(ATS_1  A)"1  ATS'1  d(  flf)  , 

~ ~ “flf  ’ 

where  x'  is  the  minimum  variance  parameter  estimate,  and 

T -1  -1 

_ (ASA) 

-x  - ~ “flf- 


(8) 

(9) 


by  Theorem  2 and  Corollary  1 in  the  section  on  Multivariate  Regression 
Analysis . 

the  correction  to  the  zeroth  estimate  is  large,  then  the 
above  procedure  should  be  iterated  because  of  the  introduction  of  the 
Taylor  series  approximation,  equation  (6). 

Since  the  radius  of  the  earth  is  sufficiently  well  known  for 
our  purposes,  we  can  restrict  our  problem  to  solving  for  position  and 
velocity  on  the  surface  of  the  assumed  earth  spheroid.  Furthermore, 
we  shall  use  a tangent  plane  system  in  place  of  the  local  spheroidal 
coordinate  system  so  that  the  study  of  the  errors  is  uniform  in  all 
directions.  Thus  the  number  of  parameters  is  reduced  to  x,  y,  x,  y 
in  the  tangent  plane . Let  x and  y designate  the  axes  in  the  tangent 
plane  where  x is  oriented  east  and  y is  directed  north.  The  data  are 
Doppler  signal  observations  which  we  assume  are  uncorrelated  random 
variables  forming  a normal  distribution  with  variance  satisfying  the 
relation 


s is  generally  unknown,  but  can  be  estimated,  and  the  weighting  fac- 
tors Pi  are  known.  Since  all  the  Doppler  signals  are  transmitted  from 
the  same  satellite  and  received  on  just  one  receiver,  we  can  set  all 
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equal  to  1.  Thus  the  weighting  matrix  R is  given  by 


where  I is  an  nxn  identity  matrix. 


Geometric  Optimization  Results 

In  order  to  obtain  geometric  optimization  results  for  satel- 
lite Doppler  navigation,  we  need  first  the  covariance  matrix  of  the 
parameters,  which  was  given  in  the  last  section  by 


S A = (aV1  A) 
“X  flf  “ 


-1 


for  a minimum  variance  estimate.  We  saw  also  that  the  weighting  matrix 
R is  given  by 


-_-o  - I > 

flf  s 


so  that  the  covariance  matrix  is  given  by 

-1 


Sa  = (AT  I A)' 

*“ S#-  r.  ^ 


“5c  s' 

- 8~(A1-A) 

Let  us  define  the  matrix  A as  follows, 


A 

nxq 


where 


ai2' 


(10) 


• * * f 


Then 
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A 

n 

— 

This  shows  that  the  covariance  matrix  of  the  parameters  is  formed  by- 
summing  sequentially  the  calculations  of  the  partial  derivatives  (the 
' 6 which  are  evaluated  at  the  zeroth  estimate). 

In  the  section  on  Probability-  Ellipsoid  of  Error  in  Chapter 
III,  we  saw  that  the  inverse  of  the  covariance  matrix,  given  in  equa- 
tion (10)  above,  is  the  equation  of  a q-dimensional  error  ellipsoid. 

If  we  take  the  upper  left  3*3  submatrix  we  obtain  the  equation  for  the 
ellipsoid  corresponding  to  error  in  position.  This  error  ellipsoid 
can  be  considered  as  centered  at  the  receiver.  If  the  equation  of  the 
ellipsoid  is  expressed  in  the  tangent  plane  system  then  the  trace  of 
the  ellipsoid  in  the  xy  - plane  is  an  ellipse  centered  at  the  receiver. 
We  obtain  different  sized  and  oriented  ellipsoids  as  we  vary  the  follow- 
ing quantities: 


i)  Standard  deviation  of  the  noise  associated  with  the 
Doppler  signal  flf; 

ii)  Number  of  observations  in  a single  satellite  pass; 
iii)  Maximum  elevation  angle  B attained  by  the  satel- 
lite with  respect  to  the  horizon  of  the  ship. 

For  different  maximum  elevation  angles  of  the  satellite  the 
error  ellipsoids  rotate  out  of  the  tangent  plane  about  the  semi -minor 
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axis.  The  relation  between  the  maximum  elevation  angle  0 to  the  sat- 
ellite and  the  elevation  angle  (J>  of  the  semi -major  axis  of  the  ellipr 
soid  for  a particular  satellite  pass,  is  given  approximately  by 

~ 700  - (7/ 9 >©-,  lo° 5 -0290° , 0c«(J)fE70o, 

where  0 and  fl)  are  in  degrees  measured  above  the  tangent  plane  and 
in  the  opposite  sense  to  each  other.  (This  relation  is  not  necessari- 
ly valid  for  receivers  at  other  than  positions  near  the  equator,  nor 
for  satellites  other  than  the  simulated  polar  satellite  given  in  the 
appendix.)  The  semi -minor  axis  remains  in  the  tangent  plane  and  is 
parallel  to  the  tangent  line  of  the  subtrack  of  the  satellite  at  the 
position  of  the  receiver.  For  the  simulated  polar  satellite  used  in 
this  analysis  the  tangent  line  to  the  subtrack  has  an  inclination  of 
about  4"  east  of  the  y-axis  of  the  tangent  plane.  This  condition  arises 
because  of  the  rotation  of  the  earth  during  the  satellite  pass  (4°  in 
16  minutes).  The  subtrack  is  defined  as  the  geoaentric  projection  of 
the  satellite  onto  the  assumed  earth  spheroid. 

The  following  table  lists  the  sizes  and  orientations  of  the 
semi -axes  of  the  error  ellipsoids,  and  the  sizes  of  those  of  the  error 
empses.>  fo**  a satellite  at  three  different  maximum  elevation  angles 
0.  The  ellipsoid  is  the  82.8  $ probability  ellipsoid  and  the  corre- 
sponding elliptical  trace  is  the  91.8  $ probability  ellipse  of  error. 
The  standard  deviation  in  the  Doppler  frequency  used  in  equation  (10) 
is 

s . - .01  cycle/sec 

ht 


relative  to  a CW  transmitted  frequency  of  f =108  mega cycles /sec . The 
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approximate  height  of  the  satellite  above  the  surface  of  the  earth 
near  the  equator  is  475  nautical  miles. 


0- 

D 

A 

C 

f 

A 

o 

O 

o\ 

o 

0 

888  km 

296.0  m 

1.4  m 

296.0  m 

-P- 

00 

0 

33° 

1139  km 

49.9  m 

1.6  m 

4.2  m 

10° 

6i 

2539  km 

222.5  m 

2.4  m 

5 .8  m 

where  -0  — Maximum  elevation  angle  of  satellite. 


<j)  — Orientation  of  major  axis  of  ellipsoid, 

D — Minimum  distance  of  satellite  from  receiver,  in  kilometers, 
A — Length  of  semi-major  axis  of  ellipsoid,  in  meters, 

C d Length  of  semi-minor  axis  of  both  ellipsoid  and  ellipse, 
in  meters, 

A Length  of  semi-major  axis  of  ellipse,  in  meters. 

We  note  that  -0-  — 90°  indicates  an  overhead  pass.  Data  are  obtained  at 
the  rate  of  one  observation  every  four  seconds,  for  a total  of  about 
220  observations  on  high  elevation  passes  and  about  200  for  low  ele- 
vation passes.  All  calculations  are  performed  on  the  IBM  7094  com- 
puter . 

The  following  figure  illustrates  the  dependence  of  the  semi- 
major axis  length  of  the  elliptical  trace  on  the  maximum  elevation 
angle  0,  attained  by  the  satellite  with  respect  to  the  receiver.  The 
variation  of  the  semi -minor  axis  length  is  much  less  and  can  be  ig- 
nored as  a measure  for  the  optimum  satellite  pass.  Thus,  if  we  define 
the  optimum  satellite  pass  as  that  pass  which  minimizes  the  semi-major 
axis  length  of  the  ellipse,  then  from  the  figure  the  optimum  satellite 
pass  appears  to  be  the  one  with  maximum  elevation  angle  0 ~ 35° • 
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< 


90°  6o°  30° 

e (degrees) 


Semi -major  axis  length  A,  in  meters,  of  ellipse  versus  maximum  ele- 
vation angle  6-,  in  degrees,  of  a satellite  pass  with  respect  to  a 
receiver  at  0°  latitude.  The  standard  deviation  of  the  Doppler  fre- 
quency s — *91  cycles  per  sec.  The  data  rate  is  one  measurement 
every  4 seconds . 
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Functional  Dependence 

Functional  dependence  among  the  errors  in  the  parameters  of 
an  estimation  problem  is  the  main  reason  for  the  need  and  development 
of  analytic  optimization  techniques . Functional  dependence  is  studied 
by  simulating  an  orbit  and,  knowing  the  correct  position  and  velocity 
of  the  receiver,  we  make  a known  error  in  the  correct  position  and 
solve  for  the  compensating  velocity  error.  Thus,  at  0°  latitude,  4o° 
maximum  elevation  angle  pass  of  the  satellite  west  of  the  receiver, 
ve  make  a known  position  error  of  Ax  = 70.7  meters,  Ay  = 70. 7 meters, 
causing  a compensating  velocity  solution  of  A x = - 57.9  cm/sec  , and 
Ay— 11.1  cm/sec. 

These  errors  are  more  easily  understood  by  studying  the  graphs1^ 
on  the  following  pages,  where  the  above  values  are  introduced  indi- 
vidually and  in  various  combinations.  In  the  graphs,  the  vertical 
measurement  is  time  (about  15  minutes),  the  horizontal  measurement  is 
the  frequency,  1 cycle/sec  full  scale.  The  vertical  line  intersecting 
the  mid -point  of  the  horizontal  axis  represents  the  Doppler  signal 
with  no  errors  in  the  parameters . Thus  these  graphs  represent  the 
errors  in  the  Doppler  signal  due  to  the  given  errors  in  the  parameters. 

The  most  striking  evidence  of  compensating  errors  is  to  note 
that  the  Doppler  error  curve  for  Ax  =70.7  meters  matches  remarkably  in 
magnitude  and  in  opposite  sign  that  for  Ay^ll.l  cm/sec..  Similarly, 
the  Doppler  error  curve  for  Ay =70.7  meters  is  equal  in  magnitude 
* 50  cm/sec  is  roughly  equal  to  1 knot  (nautical  mile/hour). 
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For  a 0 latitude,  40°  elevation  pass,  satellite  vest,  vith  combined 
errors  of  x = 70.711  m,  y = 70. 711  m,  the  velocity  solution  of  x = ~57.855 
cm/sec  and  y =11.059  cm/sec  is  obtained.  These  values  are  introduced 
individually  and  in  various  combinations  in  graphs  on  the  previous  page. 
All  four  errors  are  introduced  in  the  above  graph. 
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and  opposite  in  sign  to  that  for  Axr-  57.9  cm/sec.  When  Ax  and  Ay 
errors  are  combined,  the  resulting  error  curve  is  very  close  to  the 
true  Doppler  signal.  The  same  is  true  for  the  combination  of  fly  and 
Ax  errors,  and  for  the  final  graph  where  all  four  errors  are  combined. 

Although  the  compensating  velocity  errors  in  the  above  case 
are  not  drastically  in  error,  the  velocity  errors  can  become  as  high 
as  800  cm/sec  for  different  latitudes  and  satellite  passes. 

An  analytic  form  for  the  errors  for  a 4 0°  west  elevation  pass 
at  0 latitude  was  found  by  solving  for  compensating  velocity  errors 
for  a set  of  position  errors  of  magnitude  100  meters  at  different  ori- 
entations around  the  true  position.  An  harmonic  analysis  was  made  to 
find  the  coefficients  (kt  and  b^  ) in  the  expansion 

A x = Z k,  sin(i  O - b. ) . 
i i i 

The  values  for  ki  for  the  first  five  orders  of  Ax  are 

kx  r 77.105^ 

k2  - .0050 

kn  = .0020 

k?  — .0008 

kj.  = .0021 

A similar  result  was  obtained  for  fly.  Retaining  only  the  first  or- 
ders, therefore,  and  making  suitable  substitutions,  where 


•©■-tan  1 ( fly/  Ax), 
k =.771054  R 


where  R=100  meters. 


leads  to  the  following  form  for  Ax  and  Ay  in  terms  of  Ax  and  Ay: 

Ax  = - .048784  Ax  — .76955  Ay 
Ay  = .115695  Ax  -I-  .01*0876 Ay 

where  fix.  Ay  are  in  cm/sec  and  Ax,  Ay  in  meters. 
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Uncorrelated  Estimation  Results 


To  show  the  application  of  the  Theorem  of  Uncorrelated  Esti- 
mation we  shall  suppose  that  we  have  obtained  uncorrelated  estimates 
for  x and  y and  shall  show  the  derivation  of  equations  for  the  third 
parameters  yQ,  and  then  the  last  parameter  xQ. 

We  are  given  the  following  set  of  initial  parameters: 


With  these  values  fixed,  we  now  make  an  ordinary  minimum  variance 
estimate  for 


The  following  uncorrelated  estimates  have  been  obtained: 


*4  = y 


l t 


and  we  obtain  their  associated  covariance  matrix 


S 


x 


Thus  the  value  of  m as  given  in  the  theorem  is 


We  now  form  the  equation 


X 


2 


2 


9^ 


We  substitute  the  right-hand  side  of  equation  (ll)  for  (x| ' , x ” ) in 
the  multivariate  regression  analysis,  and  minimize  with  respect  to 
the  new  variables  y| ' , y£ ' , which  are  uncorrelated.  However,  we  must 
linearize  Af^  by  a Taylor  series  expansion  about  an  initial  estimate 
yl  > 


Thus 

A f. 


A -f*®  A f 4 . 

0fi  + i (y{' 

6 y{' 


- ^°>  + 


bflf 


6y’’ 


i (y"_  v'  '°) 
'y2  y2  ' + 


where 


Xi  - Xl)(x.  - x1)  + (Yi  - yi)(Ti  - y1)4Z1Zj 
_ ""[(Xi  - xi)2+(Yi  - yi)2+Zi2  'i 


and  where 


xi=xi’’*‘x3  - *0)  = w-  ”i’  y£’)+*$  (tt  - to)  , 

yi=xi'+fk  (h  - V = Jr2,+  Xi(  (*i  - 10)  • 


Evaluate  Af°  and  each  partial  derivative  at  y”°,  y^’°  which  are 
obtained  from 


fy’,0~ 

L 

1 

mi ' 

x°  1 

\ 

v’  ,0 
y2 

f — 

0 

1 ^ 

L-Si 

We  then  obtain  a minimum  variance  estimate  for  the  vector 
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but  we  retain  only  the  uncorrelated  estimate  y^ ' — x^'—yQ  and  fix 

its  value.  With  three  values  now  fixed,  we  solve  for  x'''  — x in 

1 o 

a one  - dimensional  regression  analysis.  The  final  uncorrelated 
estimate  solution  vector  is  then  given  by 


*}• 


The  results  obtained  here  from  the  method  of  Uncorrelated  Esti- 
mation and  for  all  the  other  methods  in  subsequent  sections  below  use 
a simulated  satellite  orbit  that  was  generated  from  equations  in  Ap- 
pendix B.  The  Doppler  signal  data  are  generated  from  equation  (4),  and 
become  random  variables  through  application  of  the  random  number  gener- 
ator in  Appendix  A,  using  a standard  deviation  in  Doppler  frequency 
sflf  =-°5  cycle/sec  and  an  initial  random  number  r^l.  The  satellite 
is  a 1+0  maximum  elevation  angle  pass  in  a polar  north  to  south  orbit 
and  is  west  of  a receiver  located  at  0°  latitude.  The  origin  of  the 
coordinate  system  is  at  the  true  position  of  the  receiver  and  the  ve- 
locity of  the  receiver  is  set  equal  to  zero  since  the  error  in  velocity 
is  independent  of  the  actual  velocity.  Thus  the  actual  position  and 
velocity  of  the  receiver  for  which  we  are  solving  in  the  multivariate 


regression  analysis  is  xQ_  0,  yQ—  0,  x_0,  y — 0.  The  regression  analy- 
sis yields  the  stationary  point  xQ=  -7.65  meters,  y0=  60.29  meters, 
x =-1+9.21  cm/sec,  y =1.63  cm/sec  , where  we  note  that  the  error  in 
velocity  is  less  than  1 knot  and  the  error  in  position  is  less  than 
the  length  of  most  ships.  This  stationary  point  is  different  for  differ- 


ent latitudes  and  maximum  elevation  angles . 
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The  method  of  Uncorrelated  Estimation  if  not  combined  with  a 
priori  constraints  consistently  yields  the  stationary  point.  Thus  we 
combine  the  two  methods  by  using  the  a priori  constraint  only  on  the 
new  variables.  This  does  not  overlap  then  with  the  results  obtained 


in  the  section  on  Constrained  Estimation.  The  position  variables  are 


in  meters  and  the  velocity  variables  are  in  cm/sec. 


Initial 
, parameter 
A \value  s 

priori  \ 

constraint  s \ 

o o 
X 

30 

70 

**  3.6 

100 

30 

70 

70.3 

100 

x° 

15 

35 

-39.2 

50 

y° 

15 

35 

11.6 

50 

Syl  Sy2 

Byj  syu 

100 

(3rd  Itei 

200 

'at ion ) 

xo 

y0 

X 

• 

y 

21.0 

-3.3 

-2.1 

1-7 

41.2 

-20.2 

11-3 

1.7 

-0.4 

51.6 

-42.5 

1.6 

56.3 
-32.8 

21.3 
1.7 

100 

(6th  Itei 

200 

'at  ion) 

xo 

X 

y 

13-1 
11.5 
-13.9 
1 .6 

23.2 

-2.3 

-3*8 

1 .6 

-3-6 

53.4 

-44.2 

1.6 

30.8 

-12.7 

3-7 

1.6 

* 

(6th  Itei 

'at  ion ) 



xo 

y0 

X 

• 

y 

22.1 

-4.7 

-2.1 

1.6 

21.1 

-16.2 

6.3 

1.6 

» The  following  rather  complicated  system  of  constraints  was 


Iteration 

1 

2 

3 

4 

5 

5 

syr  % 

100 

100 

100 

100 

100 

100 

200~ 

200 

~200| 

~200~ 

200 

200 

% 

600 

500 

4oo 

300 

200 

200 

used: 


**  Stationary  point  plus  10  units. 
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The  table  belov  is  an  example  of  the  changes  that  occur  from 
one  iteration  to  the  next.  The  solution  from  the  kth  iteration  is 
used  as  the  initial  parameter  values  for  the  k+lth  iteration.  All 
positions  are  in  meters  and  velocity  in  cm/sec  in  both  tables  below. 


Initial  parameters  x°  - y°  -'JO  m,  x°=  y°  = 35  cm/sec. 


S — s 
yl  y2 

8 H S 

y3  y4 

1 * 

2 

3 

4 

5 

6 

I 

100 

200 

*o 

yo 

X 

y 

62.9 

-6.0 

24.5 

5-9 

50.8 

-21.3 

17.4 

2.2 

41.2 
-20.2 

11.3 
1.7 

33.7 
-14.7 
5.8 
l .6 

27.9 

-8.4 

0.7 

1.6 

23.2 

-2.3 

-3.8 

1.6 

* The  numbers  1 through  6 refer  to  the  iteration  number. 


The  following  is  an  example  of  combining  Uncorrelated  and 
Constrained  Estimation  techniques  by  including  a priori  constraints 
on  both  the  new  variables  and  the  old  variables  for  all  iterations. 
The  initial  parameter  values  are  the  stationary  point  plus  100  units 
added  to  each  value . 


s 

position 

s 

velocity 

First 

Uncorrelated 

Estimate 

Second 

Uncorrelated 

Estimate 

Third 

Constrained 

Estimate 

70.6 

17.0 

6.0 

4oo 

200 

yc 

5.1 

-0.9 

2.4 

X 

-7.2 

-3-5 

-5.0 

y 

9-0 

1.8 

0.9 
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Quadratic  Pseudoinverse  Estimation  Results 
We  shall  now  apply  the  theorem  on  quadratic  pseudoinverse 
estimation  to  the  navigation  problem.  The  quadratic  minimizing  func- 
tion is  given  by 

S =||  d(  flf ) - A dx  ||jj 

where  we  can  set  the  weighting  matrix  N=Q./b^!I  and  where  the  n*4 
matrix  A will  be  assumed  to  be  of  rank  2.  This  is  a reasonable  assump- 
tion based  on  the  discussion  in  the  section  titled  Functional  Dependence. 
Thus  we  can  factor  A into  two  matrices,  each  of  rank  2,  such  that 


A = B C (12) 

nx4  nx2  2x4 


where  the  columns  of  B form  a basis  for  the  column  space  generated  by 
the  column  vectors  of  A.  Let  us  choose  as  the  basis  for  the  column 
space  the  first  two  column  vectors  of  A corresponding  to  the  position 
parameters  x and  y.  Thus  equation  (12)  becomes 


all 

a12 

a13 

al4 

an 

a12 

a2l 

a22 

a23 

a2U 

a21 

a22 

• 

• 

• 

• 

3 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

anl 

an2 

a -> 

n3 

anU 

a i 

nl 

an2 

— 

— 



The  third  column  of  A is  then  given  by 


C11  c12 
C21  c22 
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or 


(13) 


We  must  solve  equation  (13 ) for  the  vector  < CH)by  taking  the  quadratic 

( c21l 

pseudoinverse  of  B with  respect  to(]yfc%:. 


In  a similar  fashion  we  obtain 


Nj  = (|Fb)"*BT 


Sl4 

a2k 


ank 


Thus  we  have  determined  all  the  values  for 


100 


c = 


1 

0 


0 

1 


'11 


c12 
c. 


'21  "22 

Next,  define  the  inverse  of  the  covariance  matrix  for  the  parameters 

t>y 


Q = 


l/s? 


0 


1/e 


l/s2 


1/s2 * 


where  Q is  the  a priori  estimate  of  the  variances . Then  the  param- 
~T  U ^ ) 

eter  estimate  x — <xQ,  y , x,  y > is  given  by 


/ 


\ 


= 


) = Q -V  (C  ?-VT1(|‘B)-¥y  , 


lJNrl/ J^v-lJ? 


❖ 

y 

V ) 


-1  T 

where  C Q C 


2 2 2,2  2 

Sx  + C11  Sx  + c12  Sy 


2 2 

cn . c-,  s . 4 c,  ~c  s . 

11  21  x ^ 12  22  y 


cllc21  s£  + c12c22  sy 


s 4 c21  s . 4 c00  s_._ 


22 


-1  T 

and  where  Q C - 


0 

2 

c,  , s. 
11  x 


C12S, 


2 

C21sx 

c22s| 


The  following  table  contains  the  solutions  for  xQ,  y , x,  y , 
after  one  iteration  only,  since  subsequent  changes  are  negligible. 
Units  of  position  are  in  meters,  velocity  in  cm/sec. 


2 2 
sx  + a2 


x 


s 


2 

y 


a!  = -.048784,  a2=  -.76955,  b±=  .115695,  bg  — .040876  . 
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Constrained  Estimation  Results 

The  minimum  variance  constrained  estimation  for  the  parameter 
vector  x,T-(JX^  y^,  x',  y'j  is  given  by 


dx  = (FTS  F)"1  FTS  dz 


-1  ™T, 


where 


x ' - x° 

o o 


y ' y° 

d4  = x'-  x°-  l ° ° 

\ • . • n 


X' 

y' 

V 

where 
d(Af )- A f - Af°  = ) 


}*  dz  = \ 


y° 


- A f°  ' 

flfg  - Af° 


dx  ~ x 


Af  -*Af° 

' n n 


(lU) 


/ _ 


n 


i o 
xo  - xo 

^ -^o  } 

X*  - X°  ( 

r - y° 


where 


v y0*  k>  $ 


is  the  a priori  estimate  of  the  param- 
o _.,o 


T r 

eters  and  x°  — |x°,  y°,  xu,  yu  is  the  zeroth  estimate.  These  two 
vectors  often  coinside  for  the  first  solution,  given  by  x1T=)x1,  y1 

•1  -ll  * °' 

x > y j > which  is  then  used  in  the  next  regression  analysis  in  place 
of  the  zeroth  estimate . 

The  following  defines  the  F and  S matrices; 


Fr 


L 


A 

n*k 

I 

k*k 


where 


A - 


all 


21 


nl 


a12  al3  ai4 

a22  a23  a2k 


*n2 


an3  anU 
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and  S~ 


R 0 
nxn 

0 Q 
kxk 


, vhere  Rrl  I and  Q-1= 
s2 


° s£ 


si 

y 


The  matrix  Q is  the  a priori  estimate  of  the  variances  of  the 
eter . 


para  m - 


Substitute  F,  S and  dz  into  (l4)  to  give 

dx  = ( I ?AT  A + Q)_1(  i AT  d(flf ) + Q d x)  . 
& 


(15) 


y,  related 


The  following  table  lists  the  values  of  xQ,  yQ,  x, 
to  initial  parameter  values  and  a priori  constraints.  All  positions 


are  in  meters  and  velocity  in  cm/sec. 


Initial 
parameter 
values 
A \ 

priori 

constraint&\ 

x° 

P 

30 

70 

. * 

3.6 

°=  35  • 

y° 

J0 

30 

70 

70.3 

k° 

15 

35 

-39-2 

y° 

15 

35 

11.6 

sx0~  sye 

sx=sy 

6oo 

(2nd  Ite 

300 

ration ) 

xo 

y0 

X 

y 

-3.0 

4.4 

-6.3 

-0.1 

4.1 

-0.1 

-3.2 

0.5 

-5-5 

55.6 

-45*7 

1.7 

500 

(3rd  Ite 

250 

ration) 

xo 

y0 

X 

y 

-4.8 

5*1 

-6.7 

-0.3 

0.4 

0.7 

-3.6 

0.1 

-6.1 

55.6 

-45-7 

1.6 

4oo 

(4th  Ite 

200 

ration ) 

xo 

yG 

X 

y 

-4.3 

3.6 

-5.6 

-0.3 

1.4 

-0.9 

-2.5 

0.2 

-5.9 

55-5 

-45.6 

1.6 

300 

(5th  Ite 

150 

ration ) 

xo 

y0 

X 

y 

-1.2 

0.7 

-3.6 

-0.1 

8.0 

-4.2 

-0.3 

0.8 

-4.8 

55-2 

-45.4 

1.7 

Initial  parameter  values:  x°  : 

-Yn  -1 

'0;  x°~y 

S = S,r 
Xo  Vo 

s . — s . 
x y 

]_** 

2 

3 

4 

5 

6 

6oo 

300 

xo 

X 

• 

y 

25.0 

-6.2 

1.3 

2.8 

4.1 

-0.1 

-3.2 

0.5 

-4.3 

5-4 

-7.1  - 
-0.2 

-7.5 

10.2 

10.6 

-0.4 

-8.8 

14.5 

-13.8 

-0.4 

-9.2 

18.4 

-16.8 

-0.3 

* These  values  are  the  stationary  point  plus  10  units  added  to  each 
value . 

**  The  numbers  1 through  6 refer  to  the  iteration  numbers. 
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Functionally  Dependent  Constrained  Estimation  Results 

From  the  results  in  the  section  on  functional  dependence  that 
the  errors  among  the  position  and  velocity  parameters  were  of  the  form 
Ax  ~ aAx  -(-  bAy  , 

Ay  — cAx  + dAy  . 

Suppose  we  define  the  vector  of  errors  by 

Ax  = x - x , (17) 

where  x _ \xo>  yQ>  y | is  the  a priori  vector  of  parameter  esti- 


(16) 


mates 


T i • • 

» x -<>x  > Yn>  x,  y 


y (is  the  vector  of  true  values.  In  component 


form  equation  (17)  becomes 


r 

' -S 

Ax 

X - X 

0 

0 0 

Ay 
J 0 

) = ( 

- yQ 

Ax 

X - X 

h 

^ - y 

\ J 

^ J 

If  we  let  G be  defined  by 

10  0 0 


G = 


0 10  0 
a b 0 0 

c d 0 0 

where  G is  easily  seen  to  be  singular,  then  it  will  satisfy  the  relation 

x - x = GAx  , 


io6 


\ 

— 

/ 

xo 

- xo 

l 

0 

0 

0 

4xo 

- yQ 

> = 

0 

1 

0 

0 

( 

“^0 

X 

- X 

a 

b 

0 

0 

a* 

¥ 

- y 

c 

d 

0 

0 

ay 

l / 

or 


, y0  - y0 

X 


= 


a- 


a fl  xQ  + b fl  yQ 


- x 

y - y 

The  a priori  covariance  matrix  of  the  parameter  is  given  by 


c A xQ  + d A yQ 


i GS  0 
“ “fix  " 


S_  = 
x 


10  0 0 

0 10  0 

a b 0 0 

c d 0 0 


0 


0 


1 0 a c 
0 1 b d 
0 0 0 0 
0 0 0 0 


2 

2 

0 

2 

2 

sx 

sxy 

sxx 

sxy 

4 

as 

X 

csx 

s „ 

2 

s 

s . 

s • 

0 

sv 

bs2 

ds2 

xy 

y 

xy 

yy 

y 

y 

y 

s . 

s . 

s? 

s • • 

as^ 

bs^ 

a^^+b2^2 

acs2  + bds2 

XX 

yx 

X 

xy 

X 

y 

x J 

x y 

s • 

s . 

s. . 

s? 

cs2 

ds2 

2 2 
acs^  + bdSy 

22  i ,2  2 

C 8 iQ.  Sy 

X J 

xy 

yy 

xy 

y 

X 

y 

Since 

s 

X 

is  the  a 

priori 

covariance  matrix  of  the 

parameters,  we 

vert  the  above  matrix  and  substitute  it  for  Q into  equation  (15)  of  the 


previous  section  on  Constrained  Estimation  Results: 


dx  r 


4r  ATA  -I-  S 


-1 


\s 


-1 


x 


AT  d(flf  ) + S dx 
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The  main  effect  of  the  Functionally  Dependent  Constrained  Esti- 
mation method  vas  to  stabilize  the  initial  parameter  position  values 
and  to  solve  for  the  compensating  velocity  values.  Thus  a typical  re- 
sult was  the  following: 

Initial  Parameters  Results  of  Iteration  1 

.0  _ 


xo  ~ 30 


yo~  3° 
x°—  15 
Y°-  15 


xQ=  28.9 
y0  — 30.6 

x — -28.2 

y = k.S 

Subsequent  iterations  produced  negligible  changes.  Also,  the  results 

were  in  general  independent  of  the  values  of  the  a priori  constraints 

for  sx  and  s . The  only  exceptions  to  this  occurred  for  s =s  =150, 

and  sx  = l4o,  s =150.  With  initial  parameters  the  same  as  above,  the 
*/ 

results  after  6 iterations  for  sx=l4o,  sy=  150,  were 

xQ=  -12.0 

yo  - 33.2 

x = -28.2 
y = 0.02  . 
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<J)  - Plane  Navigation  Results 

One  other  method  was  developed,  but  it  is  restricted  in  its 
applicability  to  such  problems  as  two-dimensional  navigation.  Further- 
more, the  method  is  justified  only  empirically.  The  method  is  related 
to  the  plane  established  by  the  semi-minor  and  semi-major  axes  of  the 
error  ellipsoids,  where  the  elevation  angle  of  the  semi-major  axis  is 
designated  (J)  . Hence  this  plane  is  called  the  (I)  - plane.  A navi- 
gational fix  for  position  and  velocity  is  made  in  the  (J)  - plane.  The 
values  are  fixed  for  position  and  velocity  on  the  semi -minor  axis  about 
which  the  tangent  plane  has  been  rotated.  Then  a final  navigational 
fix  is  made  for  the  remaining  two  parameters  in  the  original  tangent 
plane.  The  semi -minor  axis  is  parallel  to  the  tangent  line  of  the  sub- 
track of  the  satellite  at  the  receiver's  position,  and  the  tangent  line 
is  about  4°  east  of  the  positive  y-axis,  which  is  directed  north.  In 
practice,  the  tangent  plane  is  rotated  about  the  y-axis.  Small  cor- 
rections are  made  for  the  (J)  angle  in  the  multivariate  regression  analy- 
sis by  including  a change  in  the  angle  as  an  additional  parameter  in 
the  solution.  We  note  that  the  initial  value  for  the  angle  (J)  could 
be  obtained  from  the  eigenvector  corresponding  to  the  semi -major  axis 
of  the  error  ellipsoid,  which  in  turn  can  be  obtained  directly  from 
the  raw  data.  For  more  exact  work  the  xy-axes  of  the  tangent  plane 
could  be  rotated  around  the  z-axis  (perpendicular  to  the  tangent  plane 
and  positive  when  directed  outward)  through  the  angle  corresponding  to 
the  angle  of  the  semi-minor  axis,  obtained  from  its  corresponding 
eigenve ctor . 
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In  the  table  below  are  listed  the  navigational  results  after 
six  iterations.  The  values  for  y and  y are  obtained  in  the  $ - plane. 
These  values  are  then  fixed  and  the  remaining  x and  x parameters  are 
solved  for  in  the  tangent  plane.  A priori  constraints  are  used  in 
all  cases.  All  positions  are  in  meters  and  velocity  in  cm/sec. 


%s\.  Initial 

\ parameter 

x° 

xo 

30 

70 

3.6* 

-7.64726 

A N 

^values 

y° 

•'o 

30 

70 

70.3 

60.29233 

priori  \ 

onnst.-rn  i nfX 

x° 

15 

35 

-39.2 

-49.20698 

y° 

15 

35 

11.6 

1.62542 

S ~ S 

x y 

s • — s . 

x y 

600 

300 

xo 

yQ 

X 

y 

-3.0 

-8.1 

3.2 

-0.5 

-3.3 

-9.5 

4.3 

-0.6 

-2. a 
-7.0 
2.3 

-.5 

500 

250 

xo 

X 

y 

-3.1 

-8.5 

3.5 

-0.6 

-3.4 

-10.1 

4.7 

-0.7 

-2.8 

-7.2 

2.5 

-.5 

4oo 

200 

xo 

?o 

X 

y 

-3-2 

-8.9 

3.8 

-0.6 

-3.6 

-10.7 

5.2 

-0.7 

-2.8 

-7.4 

2.6 

-0.5 

300 

150 

X 

0 

?o 

X 

y 

-3.2 

-9-2 

4.o 

-0.6 

-3.7 

-11.2 

5.6 

-0.8 

-2.8 

-7.4 

2.6 

-0.5 

150 

75 

xo 

y0 

X 

4 

y 

-0.2 

-0.8 

-2.6 

0.1 

* 


Stationary  point  plus  10  units. 
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Conclusions 

The  examples  given  in  the  tables  are  meant  to  be  representa- 
tive of  the  results  obtained  and  are  not  exhaustive.  Nevertheless,  we 
are  able  to  draw  some  general  conclusions  from  this  study  of  the  four 
derived  optimization  methods  and  the  one  empirical  method. 

The  following  methods  improved  the  estimation  results: 

Uncorrelated  Estimation 
Constrained  Estimation 
(J)  - Plane  Navigation 

The  other  two  methods,  Quadratic  Pseudoinverse  Estimation  and  Function- 
ally Dependent  Constrained  Estimation,  in  general,  did  not  improve  the 
estimate  but,  with  the  correct  choice  of  the  a priori  variance  esti- 
mate in  the  case  of  the  first  method  and  a priori  constraints  for  the 
second  method,  the  results  did  not  become  any  worse. 

The  two  methods  that  failed  are  related  to  two  others  that  the 
author  also  derived  but  are  not  presented  in  detail  here.  The  Quad- 
ratic Pseudoinverse  is  related  to  a method  that  orthogonal ized  the 
column  vectors  of  the  coefficient  matrix  of  the  parameters.  This  does 
not  improve  anything  because  the  solution  of  the  new  parameters  in  the 
transformed  space  became  the  stationary  point  when  transformed  back 
to  the  original  space.  The  Functionally  Dependent  Constrained  Esti- 
mation is  related  to  a method  that  simply  substituted  the  functional 
relations  into  the  minimizing  function,  resulting  in  a regression  analy- 
sis for  only  two  parameters.  The  answer  was  the  x,y  - components  for 
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the  stationary  point . 

Constrained  Estimation  seems  to  give  good  results  with  an  appro- 
priate choice  of  a priori  constraints  and  choice  of  iteration.  However, 
these  choices  must  be  made  before  the  regression  analysis  and  not  after- 
wards as  we  have  done , and  the  exact  dependence  must  be  determined  for 
each  problem.  We  see  also  that  if  the  initial  parameter  values  are  the 
stationary  point  or  close  to  it,  this  method  does  not  improve  the  ini- 
tial values.  We  also  see  that  further  iterations  do  not  necessarily 
improve  the  estimation. 

<P  - Plane  Navigation,  the  empirical  method,  gives  the  best  all- 
around  results . That  is,  it  is  quite  independent  of  the  choice  of  a 
priori  constraint  over  a wide  range,  and  very  neatly  handles  the 
stationary  point.  The  method  has  its  limitations,  unfortunately, 
since  it  is  limited  to  two-dimensional  space  problems  and  does  not 
seem  to  be  capable  of  generalization.  Furthermore,  since  the  ({)  - plane 
is  rotated  around  the  y-axis,  the  whole  coordinate  system  should  move 
at  the  initial  velocity  in  the  x-direction,  which  was  not  done  in  this 
example.  However,  in  one  case  which  the  author  made  with  a moving  (J)- 
plane,  the  method  produced  an  improved  estimate  and  seemed  to  have 
possibilities  of  converging. 

Uncorrelated  Estimation,  which  in  the  author's  opinion  has  the 
strongest  theoretical  justification,  did  not  produce  as  good  a set  of 
improved  values  as  the  other  two  and  also  suffered  from  the  same  prob- 
lem near  the  stationary  point  as  the  Constrained  Estimation  method. 
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However,  when  combined  with  an  a priori  constraint  on  the  old  variables 
as  well  as  the  new,  results  were  good  for  the  3rd  constrained  estimate 
after  two  uncorrelated  estimates  for  the  stationary  point  plus  100 
units  added  to  each  value . 


APPENDIX  A 


RANDOM  NUMBER  GENERATOR 


17 


The  subroutine  NDRN1  computes  a single  random  number.  A set 
or  sample  of  such  numbers  have  the  characteristics  of  a set  of  nor- 
mally distributed  numbers  with  the  given  mean  and  standard  deviation. 

To  obtain  this  set,  input  the  standard  deviation,  the  data  to 
which  random  number  is  to  be  added,  and  the  initial  random  number 
which  is  an  odd  positive,  fixed  point,  octal  integer. 

The  data  plus  random  number  is  obtained  by  the  following  for- 
mula: 


fN  =f  -f-  s (sign  (^-0.5) 


V 


ao+  al  Vi  + a2  V 


1 + b-L  Vj  + b2  V^-h  b3 


where 


fjj  _ data  plus  random  number 
f — data 

s 3 standard  deviation 

ri  - initial  random  number  (odd,  positive,  fixed  point,  octal  integer) 


Ui  — floating  point  equivalent  of  r 


Vi  = 1 - 2 lo«e  L 


aQ  = 2.515577 


a = 0.802853 
a£  = 0.010328 


0.5  (1  - 


1 - 2U1  )] 

b1  - I.U32788 
b£  - O.I89269 

b^  = 0.001308 


A different  set  of  random  numbers  can  be  obtained  by  choosing  any  new  r. . 
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APPENDIX  B 


SATELLITE  ORBIT  SIMULATION  METHOD18 


The  following  is  a description  of  the  method  used  to  gener- 
ate simulated  near  earth  satellites  which  is  considered  sufficient 
for  realistic  simulation  purposes.  The  equations  of  motion,  values 
of  the  physical  constants,  and  the  numerical  integration  program  were 
developed  by  the  Applied  Physics  Laboratory,  Johns  Hopkins  University 
(APL/JHU)  in  an  early  phase  of  their  study  of  the  Doppler  observations 
of  satellites . 

The  Force  Function  (-U,  where  U is  the  potential  function)  is 
given  by 


T 


3 6 I z 


3J, 


■+•  ■ 


22 


1- 12- 


cos  2 A 


22  | 


where  u = GM  is  the  gaussian  gravitational  constant  for  the  earth, 

(x,  y,  z)  form  a right-handed  rectangular,  inertial,  cartesian  coordi- 
nate system  where  x points  to  the  first  point  in  Aries  and  z is  the 
rotational  pole  of  the  earth;  (0,  0,  0)  lies  at  the  center  of  mass 
( = dynamical  center,  in  this  type  of  problem)  of  the  earth. 

(o)  A. 


A = arc  tan 
22  x 


'-L  y t-l 22  , 

e 57.295780  , 
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(o) 

where  A -longitude  of  Greenwich  in  radians,  measured  positive  east- 
G 

ward  from+x,  at  time  t-0,  w is  the  rotational  rate  of  the  earth,  and 


A 


22  is  a gravitational  constant  given  in  degrees  of  angular  measure . 


r “ 


2,2  2 

x + y -t-  z 


The  equations  of  motion  are : 


y=£>f^; 

oy 


Oz 


The  values  of  the  physical  constants  used  to  generate  the  ephemeris 


are 


U 

A, 

A 

A,. 


— 398,590  km3/sec2 
= 1.6237  * 10"3  x (ae)2 


3 


= 2.42X10' 


(aj: 


= 7.4  x lo“°  x (a  )k 


A 22 

A (°) 
G 

w_ 


= 1.604  xl0'5x  (a  )2 

— - 10°  or  350° 

— 1.903  038  003  radians 
= 0.7292  1159 -X  10-^  radians/second 

ae  — 6,378.166  km rr 1 earth  radius  unit  (eru) 

together  with  the  six  constants  of  integration 
xQ  - 1.1133510  eru 

yQ  = 0.31180526  eru 

z = 1.1486270  X 10  3 eru/second 


zo  — xo  — ^o  — 0. 
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The  numerical  integration  of  the  equations  of  motion  is  ac- 
complished by  a 4th  order  Runge-Kutta  program  written  by  APL/jHU.  A 
test  is  inserted  into  the  program  to  eliminate  the  ambiguity  in  longi- 
tude vhen  the  satellite  passes  over  either  pole  of  the  earth.  F F 

x'  y* 

and  Fz  are  defined  to  he  the  J22  components  of  x,  y,  and  z,  respective- 
ly. The  test,  then,  is:  If  |z/rl  = 0.99,  then  F , F and  F ~ 0;  other- 
wise  Fx,  Fy,  and  Fz  are  calculated  by  their  defined  expressions. 
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